marsh.models.new
LiaC
2018-06-05
Model selection for Ch 1 - marsh data
- Decided we should NOT include site as a random effect for all beach seine models because the low # of sites and it performed worse than regular GLMs for all.
- Early on in model selection we determined that our data has a negative binomial distribution. All models started with Poisson distributions and all were overdispersed, which was corrected by the neg binomial distribution. All models here use negative binomial distribution with a log link.
- We used AICc to rank all models in dredge, as it penalizes the strength of likelihood more for very small sample sizes, and is more adaptive to moderate sample sizes, so that we could apply the same information criterion across all models (Harrison et al 2018; Brewer et al 2016) – could explain better…
- We used model averaging following a delta AICc of less than 4
Chinook beach seine model
response = average abundance of Chinook salmon [per beach seine site]
This is 1 of 4 separate marsh models for (1)Chinook, (2)chum, (3)other migratory species and (4)resident species
Load all data - standardize variables, create subgroups from beach seine data
# load Beach seine data aggregated by species/set; note this contains
# unidentified larval fish
b <- read.csv("/Users/Lia/Documents/Git/Fraser-salmon/all.data/beach.catch.csv")
summary(b)## Year Date Temp.surf DO.surf
## Min. :0.0000 5/17/2017: 42 Min. : 3.43 Min. : 48.80
## 1st Qu.:0.0000 30-May-16: 35 1st Qu.:11.68 1st Qu.: 90.90
## Median :0.0000 5/31/2017: 32 Median :13.30 Median :103.20
## Mean :0.4817 6/16/2017: 31 Mean :13.36 Mean : 99.55
## 3rd Qu.:1.0000 23-Sep-16: 24 3rd Qu.:16.30 3rd Qu.:112.80
## Max. :1.0000 7/12/2017: 24 Max. :18.32 Max. :136.60
## (Other) :632
## DOmg.surf pH.surf Sal.surf Temp.bot
## Min. : 4.68 Min. : 6.440 Min. : 0.000 Min. : 4.15
## 1st Qu.: 9.07 1st Qu.: 7.710 1st Qu.: 0.110 1st Qu.:11.58
## Median :10.74 Median : 7.905 Median : 0.890 Median :13.20
## Mean :10.46 Mean : 8.059 Mean : 2.326 Mean :13.36
## 3rd Qu.:11.72 3rd Qu.: 8.150 3rd Qu.: 3.270 3rd Qu.:16.30
## Max. :15.50 Max. :10.460 Max. :19.400 Max. :18.32
##
## DO.bot DOmg.bot pH.bot Sal.bot
## Min. : 48.50 Min. : 4.68 Min. :6.440 Min. : 0.050
## 1st Qu.: 93.20 1st Qu.: 9.21 1st Qu.:7.750 1st Qu.: 0.120
## Median :101.40 Median :10.49 Median :7.895 Median : 0.850
## Mean : 99.13 Mean :10.40 Mean :8.021 Mean : 2.416
## 3rd Qu.:112.50 3rd Qu.:11.54 3rd Qu.:8.150 3rd Qu.: 3.340
## Max. :140.30 Max. :15.50 Max. :9.890 Max. :20.750
##
## J.date SIN_TIME COS_TIME Round
## Min. : 64.0 Min. :-0.9929 Min. :-0.9990 Min. : 1.000
## 1st Qu.:124.0 1st Qu.:-0.1436 1st Qu.:-0.9641 1st Qu.: 4.000
## Median :151.0 Median : 0.5176 Median :-0.7075 Median : 6.000
## Mean :157.2 Mean : 0.3478 Mean :-0.6395 Mean : 5.901
## 3rd Qu.:191.0 3rd Qu.: 0.8460 3rd Qu.:-0.4431 3rd Qu.: 8.000
## Max. :282.0 Max. : 1.0000 Max. : 0.4527 Max. :10.000
##
## Habitat Site Set Species
## Marsh:820 M1:155 Min. :1.000 Three-spined stickleback:161
## M2:179 1st Qu.:1.000 Chinook :124
## M3:171 Median :2.000 Peamouth chub :119
## M4:172 Mean :2.059 Starry flounder : 98
## M5:143 3rd Qu.:3.000 Staghorn sculpin : 74
## Max. :7.000 Prickly sculpin : 50
## (Other) :194
## Life_stage abundance class round
## 0 : 28 Min. : 0.00 0 : 28 Min. : 0.500
## adult :176 1st Qu.: 1.00 migratory:201 1st Qu.: 4.000
## juvenile:474 Median : 2.00 resident :588 Median : 6.000
## NA's :142 Mean : 8.69 NA's : 3 Mean : 5.757
## 3rd Qu.: 6.00 3rd Qu.: 8.000
## Max. :509.00 Max. :10.000
##
## month
## May :210
## June :163
## April :160
## July :119
## March : 81
## September: 45
## (Other) : 42sapply(b, sd)## Year Date Temp.surf DO.surf DOmg.surf pH.surf
## 0.4999702 16.0088878 3.2129387 19.5469879 2.0540975 0.6425623
## Sal.surf Temp.bot DO.bot DOmg.bot pH.bot Sal.bot
## 3.3653542 3.1011628 19.8481953 2.1137164 0.5439179 3.6501406
## J.date SIN_TIME COS_TIME Round Habitat Site
## 47.0456528 0.5773336 0.3706291 2.5508623 0.0000000 1.3720692
## Set Species Life_stage abundance class round
## 0.9564320 6.3846921 NA 29.7143774 NA 2.2935398
## month
## 2.0433562### Separate into groups; automatically removes 3 x unidentified larval and 28
### x '0' abundance rows 1
b.1 <- b[which(b$Species == "Chinook"), ]
### standardize continuous variables
library(robustHD)
# standardize variables to be centered on the mean (mean becomes 0) using
# the standardize function from robustHD
b.1$s.temp <- standardize(b.1$Temp.surf, centerFun = mean, scaleFun = sd)
b.1$s.sal <- standardize(b.1$Sal.surf, centerFun = mean, scaleFun = sd)
b.1$s.do <- standardize(b.1$DOmg.surf, centerFun = mean, scaleFun = sd)
b.1$s.pH <- standardize(b.1$pH.surf, centerFun = mean, scaleFun = sd)
b.1$s.J.date <- standardize(b.1$J.date, centerFun = mean, scaleFun = sd)
### Create variable j2 which is Julian day squared in order to represent J
### date as quadratic relationship instead of linear
b.1$j2 <- b.1$s.J.date^2
# 1 summarize by site-day
library(plyr)
b.chin <- ddply(b.1, .(Year, J.date, s.J.date, j2, Habitat, Site, s.temp, s.sal,
s.do, s.pH), summarize, abundance = sum(abundance))
# plot abundance over julian day to see distribution:
plot(abundance ~ J.date, data = b.chin)
Habitat variables
Load habitat variables – these are habitat characteristics measured at each marsh site once. See “habitat characteristics” section of model.b3.html for full descriptions
## Add marsh habitat variables for each site
hab <- read.csv("/Users/Lia/Documents/Git/Fraser-salmon/all.data/site.char.marsh.csv")
hab2 <- hab[, c(1, 6, 8:10, 13, 19, 20)]
summary(hab2)## Site stnwidth vegheight shtdensity shtdenhigh
## :13 Min. : 18.80 Min. :0.3235 Min. : 176.0 Min. : 752
## M1: 4 1st Qu.: 24.30 1st Qu.:0.6747 1st Qu.: 344.0 1st Qu.:1664
## M2: 1 Median : 53.90 Median :1.2538 Median : 709.3 Median :1837
## M3: 4 Mean : 51.26 Mean :1.2451 Mean : 812.2 Mean :1843
## M4: 4 3rd Qu.: 60.30 3rd Qu.:1.5718 3rd Qu.:1392.0 3rd Qu.:1933
## M5: 4 Max. :170.07 Max. :2.5837 Max. :1792.0 Max. :2768
## NA's :13 NA's :13 NA's :13 NA's :13
## tcelev angbank meanturb
## Min. :0.3440 Min. : 4.000 Min. :29.87
## 1st Qu.:0.6390 1st Qu.: 8.531 1st Qu.:33.67
## Median :0.8470 Median : 9.500 Median :37.47
## Mean :0.8408 Mean :12.144 Mean :35.95
## 3rd Qu.:1.0500 3rd Qu.:12.917 3rd Qu.:41.02
## Max. :1.3590 Max. :31.000 Max. :41.52
## NA's :13 NA's :13 NA's :13sapply(hab2, sd, na.rm = TRUE)## Warning in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm = na.rm): Calling var(x) on a factor x is deprecated and will become an error.
## Use something like 'all(duplicated(x)[-1L])' to test for a constant vector.## Site stnwidth vegheight shtdensity shtdenhigh tcelev
## 1.9546584 37.6011569 0.6223313 552.9042986 475.4831598 0.2944439
## angbank meanturb
## 6.6741386 4.5283513b.chin.hab <- b.chin
b.chin.hab$stnwidth <- hab2[match(b.chin.hab$Site, hab$Site), 2] #add stnwidth to b.chin.hab by matching with Site #. Then repeat with all other habitat variables
b.chin.hab$vegelev <- hab2[match(b.chin.hab$Site, hab$Site), 3]
b.chin.hab$shtdensity <- hab2[match(b.chin.hab$Site, hab$Site), 4]
b.chin.hab$shtdenhigh <- hab2[match(b.chin.hab$Site, hab$Site), 5]
b.chin.hab$tcelev <- hab2[match(b.chin.hab$Site, hab$Site), 6]
b.chin.hab$angbank <- hab2[match(b.chin.hab$Site, hab$Site), 7]
b.chin.hab$meanturb <- hab2[match(b.chin.hab$Site, hab$Site), 8]
## standardize habitat variables using 'standardize' function from robustHD
## package
b.chin.hab$stnwidth <- standardize(b.chin.hab$stnwidth, centerFun = mean, scaleFun = sd)
b.chin.hab$stnwidth <- as.numeric(b.chin.hab$stnwidth)
b.chin.hab$vegelev <- standardize(b.chin.hab$vegelev, centerFun = mean, scaleFun = sd)
b.chin.hab$vegelev <- as.numeric(b.chin.hab$vegelev)
b.chin.hab$shtdensity <- standardize(b.chin.hab$shtdensity, centerFun = mean,
scaleFun = sd)
b.chin.hab$shtdensity <- as.numeric(b.chin.hab$shtdensity)
b.chin.hab$shtdenhigh <- standardize(b.chin.hab$shtdenhigh, centerFun = mean,
scaleFun = sd)
b.chin.hab$shtdenhigh <- as.numeric(b.chin.hab$shtdenhigh)
b.chin.hab$tcelev <- standardize(b.chin.hab$tcelev, centerFun = mean, scaleFun = sd)
b.chin.hab$tcelev <- as.numeric(b.chin.hab$tcelev)
b.chin.hab$angbank <- standardize(b.chin.hab$angbank, centerFun = mean, scaleFun = sd)
b.chin.hab$angbank <- as.numeric(b.chin.hab$angbank)
b.chin.hab$meanturb <- standardize(b.chin.hab$meanturb, centerFun = mean, scaleFun = sd)
b.chin.hab$meanturb <- as.numeric(b.chin.hab$meanturb)VIF for collinearity of habitat variables
Assess variance inflation factors
Note that despite high collinearity of temp, when we tried to remove it the model would not converge. When we ran just temperature and removed julian day and j2, the AIC and deviance explained (r2) were considerably worse. We determined in this case it would be best to leave it in.
library(car)
library(GGally)
# alias function identifies covariates that are multiples of each other - in
# this case have some habitat parameters causing issues.
alias(lm(abundance ~ s.J.date + j2 + Year + s.temp + s.sal + s.do + s.pH + stnwidth +
vegelev + shtdenhigh + shtdensity + tcelev + angbank + meanturb, data = b.chin.hab))
## Model :
## abundance ~ s.J.date + j2 + Year + s.temp + s.sal + s.do + s.pH +
## stnwidth + vegelev + shtdenhigh + shtdensity + tcelev + angbank +
## meanturb
##
## Complete :
## (Intercept) s.J.date j2 Year
## tcelev 0 0 0 0
## angbank 0 0 0 0
## meanturb 0 0 0 0
## s.temp s.sal s.do s.pH
## tcelev 0 0 0 0
## angbank 0 0 0 0
## meanturb 0 0 0 0
## stnwidth vegelev shtdenhigh shtdensity
## tcelev 17695/31146 16003/35591 -19915/37564 1770469/5798310
## angbank -58882/110101 258521/352474 -1628/3495 -20031/306857
## meanturb 3877/11296 -48133/68263 -14929/23063 -1997/4121
# We determined hab variables for Chinook should be reduced to stnwidth,
# meanturb, and vegelev based on significance and collinearity when doing
# model exploration (see marsh.chinook.Rmd).
vif(lm(abundance ~ s.J.date + j2 + Year + s.temp + s.sal + s.do + s.pH + stnwidth +
vegelev + meanturb, data = b.chin.hab))
## s.J.date j2 Year s.temp s.sal s.do s.pH stnwidth
## 6.572784 1.992816 1.802451 6.767715 2.158658 2.514453 1.512845 1.920460
## vegelev meanturb
## 1.726731 1.642707
## Pearson Corr with all vars
year <- b.chin.hab$Year
jday <- b.chin.hab$s.J.date
j2 <- b.chin.hab$j2
temp <- b.chin.hab$s.temp
do <- b.chin.hab$s.do
ph <- b.chin.hab$s.pH
sal <- b.chin.hab$s.sal
veg <- b.chin.hab$vegelev
turb <- b.chin.hab$meanturb
width <- b.chin.hab$stnwidth
shootden <- b.chin.hab$shtdensity
shtdenhigh <- b.chin.hab$shtdenhigh
elev <- b.chin.hab$tcelev
angbank <- b.chin.hab$angbank
habcovar <- cbind.data.frame(year, jday, j2, temp, do, ph, sal, veg, turb, width,
shootden, shtdenhigh, elev, angbank)
ggpairs(data = na.omit(habcovar), title = "Pearson Correlation plot for all variables in beach seine Chinook salmon dataset")
ggsave("/Users/Lia/Documents/Git/Fraser-salmon/Ch1.Seasonal_diversity/expl.plots/PearsonCorrAllVars_chinookmarsh.pdf")Model selection
Full model
The beach seine models performed better without site as a random effect (e.g. see marsh.chinook.Rmd), and only had a maximum of 5 sites in each data set, so we decided not to include site as a RE in the global models.
library(MASS)
## Chinook: full model with top habitat variables
chin <- glm.nb(abundance ~ s.J.date + j2 + Year + s.temp + s.sal + s.do + s.pH +
stnwidth + vegelev + meanturb, data = b.chin.hab, na.action = "na.fail")
summary(chin) #AIC 441.56, deviance explained = 100*((132.69-64.26)/132.69) = 51.57 % ((note this is same calculation that function below uses for pseudo r-squared))##
## Call:
## glm.nb(formula = abundance ~ s.J.date + j2 + Year + s.temp +
## s.sal + s.do + s.pH + stnwidth + vegelev + meanturb, data = b.chin.hab,
## na.action = "na.fail", init.theta = 1.448876434, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4538 -0.9538 -0.4630 0.4737 1.9198
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.857958 0.231350 8.031 9.67e-16 ***
## s.J.date -1.528136 0.270307 -5.653 1.57e-08 ***
## j2 -0.074091 0.126257 -0.587 0.557320
## Year 1.043695 0.309056 3.377 0.000733 ***
## s.temp 0.927479 0.274609 3.377 0.000732 ***
## s.sal -0.097340 0.161957 -0.601 0.547825
## s.do -0.115719 0.174303 -0.664 0.506757
## s.pH 0.210129 0.146060 1.439 0.150249
## stnwidth 0.372425 0.157667 2.362 0.018172 *
## vegelev -0.044583 0.150723 -0.296 0.767389
## meanturb 0.002958 0.148951 0.020 0.984157
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.4489) family taken to be 1)
##
## Null deviance: 132.69 on 62 degrees of freedom
## Residual deviance: 64.26 on 52 degrees of freedom
## AIC: 441.56
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.449
## Std. Err.: 0.279
##
## 2 x log-likelihood: -417.559chin2 <- glm.nb(abundance ~ s.J.date + j2 + Year + s.temp + s.sal + s.do + s.pH +
shtdensity + angbank + tcelev, data = b.chin.hab, na.action = "na.fail")
summary(chin2) #AIC 441.55##
## Call:
## glm.nb(formula = abundance ~ s.J.date + j2 + Year + s.temp +
## s.sal + s.do + s.pH + shtdensity + angbank + tcelev, data = b.chin.hab,
## na.action = "na.fail", init.theta = 1.452371266, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4044 -0.9797 -0.4175 0.5125 2.0049
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.85729 0.23027 8.066 7.28e-16 ***
## s.J.date -1.51924 0.26942 -5.639 1.71e-08 ***
## j2 -0.07838 0.12493 -0.627 0.530397
## Year 1.05217 0.31102 3.383 0.000717 ***
## s.temp 0.93394 0.27886 3.349 0.000811 ***
## s.sal -0.08800 0.16035 -0.549 0.583136
## s.do -0.11516 0.17403 -0.662 0.508158
## s.pH 0.20512 0.14431 1.421 0.155207
## shtdensity -0.10413 0.12897 -0.807 0.419426
## angbank -0.33802 0.15674 -2.156 0.031046 *
## tcelev 0.37826 0.16659 2.271 0.023166 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.4524) family taken to be 1)
##
## Null deviance: 132.965 on 62 degrees of freedom
## Residual deviance: 64.384 on 52 degrees of freedom
## AIC: 441.55
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.452
## Std. Err.: 0.280
##
## 2 x log-likelihood: -417.554vif(chin2)## s.J.date j2 Year s.temp s.sal s.do
## 6.155167 1.967313 1.844052 6.557370 2.105972 2.502986
## s.pH shtdensity angbank tcelev
## 1.430439 1.251839 1.881877 2.102137## Now look at the variance inflation factors of this full model
vif(glm.nb(abundance ~ s.J.date + j2 + Year + s.temp + s.sal + s.do + s.pH +
stnwidth + vegelev + meanturb, data = b.chin.hab))## s.J.date j2 Year s.temp s.sal s.do s.pH stnwidth
## 6.182463 2.004886 1.816737 6.348249 2.145333 2.501217 1.462386 1.937416
## vegelev meanturb
## 1.706736 1.669399# Test fit and assumptions of full (global) model - Jday is high (6.18)
library(sjPlot)## Warning: package 'sjPlot' was built under R version 3.4.3## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.2.12
## Current Matrix version is 1.2.11
## Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package# Check model assumptions - note that we are not actually assuming linear
# relationships with Jday, but can't tell this general function that.
sjp.glm(chin, type = "ma")## No outliers detected.## NULL
## --------------------
## Check significance of terms when they entered the model...## Anova:## Warning in anova.negbin(logreg, test = "Chisq"): tests made without re-
## estimating 'theta'
## Analysis of Deviance Table
##
## Model: Negative Binomial(1.4489), link: log
##
## Response: abundance
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 62 132.688
## s.J.date 1 43.330 61 89.359 4.625e-11 ***
## j2 1 1.013 60 88.346 0.3142812
## Year 1 0.775 59 87.571 0.3786187
## s.temp 1 10.904 58 76.667 0.0009597 ***
## s.sal 1 0.271 57 76.397 0.6028306
## s.do 1 1.269 56 75.128 0.2599757
## s.pH 1 1.343 55 73.784 0.2464322
## stnwidth 1 9.435 54 64.349 0.0021283 **
## vegelev 1 0.089 53 64.260 0.7660493
## meanturb 1 0.000 52 64.260 0.9849889
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# See pseudo R2 for full model
library(piecewiseSEM)
rsquared(chin, aicc = TRUE)## Class Family Link n R.squared
## 1 negbin Negative Binomial log 63 0.515708Model selection - dredge
Then continued to dredge function to determine optimal model using AICc (chin has 63 observations)
library(MuMIn)
library(knitr)
# Generate model set
model.set.chin <- dredge(chin)
# Create model selection table
model_table.chin <- model.sel(model.set.chin)
options(scipen = 7)
names(model_table.chin) <- c("(Int)", "Jday^2", "Mean turb", "DO", "Jday", "pH",
"Sal", "Temp", "Width", "Veg. elev", "Year", "df", "LL", "AICc", "delta",
"weight")
kable(head(model_table.chin, n = 100), digits = 3)| (Int) | Jday^2 | Mean turb | DO | Jday | pH | Sal | Temp | Width | Veg. elev | Year | df | LL | AICc | delta | weight | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 729 | 1.839 | NA | NA | NA | -1.449 | 0.174 | NA | 0.908 | 0.323 | NA | 0.959 | 7 | -209.394 | 434.824 | 0.000 | 0.1139834028 |
| 713 | 1.834 | NA | NA | NA | -1.418 | NA | NA | 0.886 | 0.306 | NA | 1.003 | 6 | -210.904 | 435.308 | 0.484 | 0.0894893114 |
| 969 | 1.822 | NA | NA | NA | -1.441 | NA | NA | 0.905 | 0.355 | -0.130 | 1.014 | 7 | -210.387 | 436.810 | 1.986 | 0.0422265266 |
| 761 | 1.801 | NA | NA | NA | -1.471 | 0.176 | -0.120 | 0.889 | 0.335 | NA | 1.026 | 8 | -209.079 | 436.825 | 2.001 | 0.0419136723 |
| 985 | 1.831 | NA | NA | NA | -1.456 | 0.158 | NA | 0.913 | 0.348 | -0.069 | 0.971 | 8 | -209.258 | 437.183 | 2.359 | 0.0350370532 |
| 733 | 1.848 | NA | NA | -0.063 | -1.463 | 0.201 | NA | 0.895 | 0.337 | NA | 0.929 | 8 | -209.283 | 437.233 | 2.409 | 0.0341737131 |
| 745 | 1.800 | NA | NA | NA | -1.438 | NA | -0.111 | 0.869 | 0.314 | NA | 1.064 | 7 | -210.675 | 437.386 | 2.562 | 0.0316545586 |
| 715 | 1.831 | NA | 0.091 | NA | -1.396 | NA | NA | 0.866 | 0.265 | NA | 1.001 | 7 | -210.683 | 437.402 | 2.579 | 0.0313987993 |
| 731 | 1.838 | NA | 0.026 | NA | -1.442 | 0.169 | NA | 0.901 | 0.311 | NA | 0.961 | 8 | -209.377 | 437.422 | 2.598 | 0.0310983070 |
| 730 | 1.848 | -0.016 | NA | NA | -1.456 | 0.173 | NA | 0.920 | 0.323 | NA | 0.976 | 8 | -209.380 | 437.427 | 2.603 | 0.0310139791 |
| 717 | 1.827 | NA | NA | 0.062 | -1.410 | NA | NA | 0.904 | 0.296 | NA | 1.026 | 7 | -210.763 | 437.562 | 2.738 | 0.0289933656 |
| 714 | 1.853 | -0.034 | NA | NA | -1.432 | NA | NA | 0.912 | 0.306 | NA | 1.036 | 7 | -210.845 | 437.726 | 2.903 | 0.0267000550 |
| 587 | 1.884 | NA | 0.242 | NA | -1.334 | NA | NA | 0.821 | NA | NA | 0.946 | 6 | -212.644 | 438.788 | 3.964 | 0.0157038322 |
| 970 | 1.842 | -0.038 | NA | NA | -1.458 | NA | NA | 0.935 | 0.356 | -0.133 | 1.053 | 8 | -210.312 | 439.291 | 4.467 | 0.0122144611 |
| 973 | 1.818 | NA | NA | 0.036 | -1.435 | NA | NA | 0.914 | 0.346 | -0.121 | 1.028 | 8 | -210.343 | 439.352 | 4.529 | 0.0118423972 |
| 971 | 1.821 | NA | 0.041 | NA | -1.429 | NA | NA | 0.895 | 0.331 | -0.115 | 1.013 | 8 | -210.350 | 439.367 | 4.544 | 0.0117550191 |
| 1001 | 1.810 | NA | NA | NA | -1.446 | NA | -0.044 | 0.896 | 0.353 | -0.114 | 1.037 | 8 | -210.359 | 439.384 | 4.560 | 0.0116560358 |
| 765 | 1.810 | NA | NA | -0.055 | -1.482 | 0.200 | -0.116 | 0.878 | 0.347 | NA | 0.996 | 9 | -208.997 | 439.391 | 4.567 | 0.0116169357 |
| 762 | 1.813 | -0.026 | NA | NA | -1.484 | 0.175 | -0.125 | 0.909 | 0.335 | NA | 1.057 | 9 | -209.043 | 439.482 | 4.658 | 0.0111019226 |
| 763 | 1.799 | NA | 0.030 | NA | -1.463 | 0.171 | -0.121 | 0.881 | 0.321 | NA | 1.028 | 9 | -209.057 | 439.509 | 4.686 | 0.0109485748 |
| 747 | 1.794 | NA | 0.097 | NA | -1.417 | NA | -0.115 | 0.848 | 0.272 | NA | 1.065 | 8 | -210.425 | 439.516 | 4.692 | 0.0109115390 |
| 1017 | 1.803 | NA | NA | NA | -1.471 | 0.171 | -0.109 | 0.892 | 0.341 | -0.020 | 1.023 | 9 | -209.070 | 439.537 | 4.713 | 0.0108011948 |
| 749 | 1.787 | NA | NA | 0.073 | -1.431 | NA | -0.121 | 0.889 | 0.304 | NA | 1.098 | 8 | -210.481 | 439.629 | 4.805 | 0.0103155662 |
| 989 | 1.840 | NA | NA | -0.069 | -1.471 | 0.187 | NA | 0.898 | 0.365 | -0.074 | 0.938 | 9 | -209.126 | 439.648 | 4.824 | 0.0102158399 |
| 734 | 1.886 | -0.060 | NA | -0.109 | -1.499 | 0.217 | NA | 0.932 | 0.347 | NA | 0.971 | 9 | -209.149 | 439.694 | 4.871 | 0.0099816884 |
| 746 | 1.821 | -0.043 | NA | NA | -1.459 | NA | -0.118 | 0.902 | 0.316 | NA | 1.112 | 8 | -210.580 | 439.827 | 5.003 | 0.0093401344 |
| 719 | 1.825 | NA | 0.081 | 0.049 | -1.392 | NA | NA | 0.882 | 0.262 | NA | 1.021 | 8 | -210.597 | 439.861 | 5.038 | 0.0091813130 |
| 986 | 1.843 | -0.021 | NA | NA | -1.466 | 0.157 | NA | 0.929 | 0.348 | -0.072 | 0.995 | 9 | -209.234 | 439.864 | 5.040 | 0.0091695253 |
| 987 | 1.831 | NA | 0.000 | NA | -1.456 | 0.159 | NA | 0.913 | 0.348 | -0.069 | 0.971 | 9 | -209.258 | 439.913 | 5.089 | 0.0089497627 |
| 735 | 1.846 | NA | 0.026 | -0.063 | -1.456 | 0.196 | NA | 0.888 | 0.325 | NA | 0.930 | 9 | -209.267 | 439.929 | 5.106 | 0.0088745728 |
| 716 | 1.845 | -0.025 | 0.087 | NA | -1.407 | NA | NA | 0.886 | 0.267 | NA | 1.026 | 8 | -210.650 | 439.968 | 5.144 | 0.0087068344 |
| 603 | 1.892 | NA | 0.214 | NA | -1.357 | 0.127 | NA | 0.836 | NA | NA | 0.919 | 7 | -211.997 | 440.030 | 5.206 | 0.0084387393 |
| 585 | 1.937 | NA | NA | NA | -1.373 | NA | NA | 0.855 | NA | NA | 0.904 | 5 | -214.515 | 440.083 | 5.259 | 0.0082185116 |
| 732 | 1.846 | -0.015 | 0.024 | NA | -1.449 | 0.169 | NA | 0.912 | 0.311 | NA | 0.976 | 9 | -209.366 | 440.129 | 5.305 | 0.0080338508 |
| 718 | 1.832 | -0.008 | NA | 0.057 | -1.414 | NA | NA | 0.908 | 0.296 | NA | 1.032 | 8 | -210.761 | 440.188 | 5.364 | 0.0077990801 |
| 601 | 1.941 | NA | NA | NA | -1.393 | 0.163 | NA | 0.873 | NA | NA | 0.868 | 6 | -213.444 | 440.388 | 5.565 | 0.0070544988 |
| 619 | 1.853 | NA | 0.247 | NA | -1.355 | NA | -0.097 | 0.807 | NA | NA | 1.001 | 7 | -212.468 | 440.972 | 6.148 | 0.0052690565 |
| 591 | 1.875 | NA | 0.223 | 0.070 | -1.330 | NA | NA | 0.845 | NA | NA | 0.974 | 7 | -212.483 | 441.002 | 6.178 | 0.0051910483 |
| 843 | 1.883 | NA | 0.246 | NA | -1.327 | NA | NA | 0.814 | NA | 0.039 | 0.945 | 7 | -212.602 | 441.241 | 6.417 | 0.0046064972 |
| 588 | 1.892 | -0.015 | 0.240 | NA | -1.341 | NA | NA | 0.834 | NA | NA | 0.960 | 7 | -212.634 | 441.305 | 6.481 | 0.0044615718 |
| 589 | 1.916 | NA | NA | 0.129 | -1.359 | NA | NA | 0.893 | NA | NA | 0.955 | 6 | -213.965 | 441.431 | 6.607 | 0.0041892510 |
| 217 | 2.396 | NA | NA | NA | -1.017 | 0.186 | NA | 0.394 | 0.289 | NA | NA | 6 | -214.040 | 441.580 | 6.756 | 0.0038884300 |
| 766 | 1.850 | -0.068 | NA | -0.106 | -1.525 | 0.218 | -0.123 | 0.919 | 0.358 | NA | 1.049 | 10 | -208.824 | 441.879 | 7.055 | 0.0033489426 |
| 1002 | 1.830 | -0.042 | NA | NA | -1.467 | NA | -0.053 | 0.928 | 0.353 | -0.113 | 1.085 | 9 | -210.270 | 441.936 | 7.112 | 0.0032541919 |
| 1005 | 1.799 | NA | NA | 0.047 | -1.440 | NA | -0.062 | 0.904 | 0.339 | -0.095 | 1.065 | 9 | -210.289 | 441.974 | 7.150 | 0.0031927535 |
| 972 | 1.840 | -0.035 | 0.032 | NA | -1.447 | NA | NA | 0.924 | 0.337 | -0.121 | 1.048 | 9 | -210.290 | 441.976 | 7.153 | 0.0031890868 |
| 1003 | 1.803 | NA | 0.058 | NA | -1.432 | NA | -0.066 | 0.877 | 0.317 | -0.085 | 1.046 | 9 | -210.292 | 441.980 | 7.156 | 0.0031834838 |
| 751 | 1.784 | NA | 0.084 | 0.059 | -1.414 | NA | -0.122 | 0.868 | 0.268 | NA | 1.094 | 9 | -210.300 | 441.996 | 7.172 | 0.0031576250 |
| 974 | 1.837 | -0.032 | NA | 0.014 | -1.453 | NA | NA | 0.934 | 0.352 | -0.129 | 1.052 | 9 | -210.308 | 442.011 | 7.188 | 0.0031338963 |
| 975 | 1.817 | NA | 0.037 | 0.033 | -1.424 | NA | NA | 0.904 | 0.325 | -0.108 | 1.026 | 9 | -210.313 | 442.023 | 7.199 | 0.0031158397 |
| 990 | 1.886 | -0.075 | NA | -0.128 | -1.518 | 0.204 | NA | 0.945 | 0.383 | -0.089 | 0.992 | 10 | -208.927 | 442.084 | 7.260 | 0.0030217716 |
| 748 | 1.811 | -0.034 | 0.090 | NA | -1.435 | NA | -0.120 | 0.875 | 0.275 | NA | 1.102 | 9 | -210.369 | 442.133 | 7.310 | 0.0029482865 |
| 767 | 1.808 | NA | 0.030 | -0.055 | -1.474 | 0.194 | -0.117 | 0.869 | 0.332 | NA | 0.999 | 10 | -208.975 | 442.180 | 7.357 | 0.0028797643 |
| 1021 | 1.813 | NA | NA | -0.058 | -1.482 | 0.194 | -0.098 | 0.881 | 0.356 | -0.030 | 0.990 | 10 | -208.978 | 442.188 | 7.364 | 0.0028693484 |
| 859 | 1.894 | NA | 0.219 | NA | -1.348 | 0.146 | NA | 0.825 | NA | 0.087 | 0.910 | 8 | -211.793 | 442.254 | 7.430 | 0.0027763221 |
| 764 | 1.811 | -0.024 | 0.026 | NA | -1.476 | 0.170 | -0.125 | 0.900 | 0.322 | NA | 1.057 | 10 | -209.025 | 442.282 | 7.458 | 0.0027378188 |
| 635 | 1.862 | NA | 0.220 | NA | -1.377 | 0.126 | -0.096 | 0.821 | NA | NA | 0.973 | 8 | -211.810 | 442.286 | 7.463 | 0.0027311620 |
| 1018 | 1.815 | -0.027 | NA | NA | -1.484 | 0.169 | -0.112 | 0.913 | 0.342 | -0.022 | 1.054 | 10 | -209.032 | 442.294 | 7.471 | 0.0027203250 |
| 1019 | 1.800 | NA | 0.027 | NA | -1.464 | 0.169 | -0.116 | 0.883 | 0.325 | -0.008 | 1.027 | 10 | -209.055 | 442.341 | 7.518 | 0.0026570316 |
| 750 | 1.796 | -0.015 | NA | 0.064 | -1.439 | NA | -0.122 | 0.898 | 0.305 | NA | 1.110 | 9 | -210.473 | 442.342 | 7.519 | 0.0026557726 |
| 617 | 1.914 | NA | NA | NA | -1.387 | NA | -0.076 | 0.842 | NA | NA | 0.946 | 6 | -214.424 | 442.347 | 7.524 | 0.0026490612 |
| 586 | 1.955 | -0.035 | NA | NA | -1.390 | NA | NA | 0.887 | NA | NA | 0.941 | 6 | -214.463 | 442.426 | 7.602 | 0.0025472267 |
| 201 | 2.421 | NA | NA | NA | -0.969 | NA | NA | 0.341 | 0.266 | NA | NA | 5 | -215.700 | 442.453 | 7.629 | 0.0025131666 |
| 991 | 1.840 | NA | -0.002 | -0.069 | -1.472 | 0.187 | NA | 0.899 | 0.367 | -0.075 | 0.938 | 10 | -209.126 | 442.482 | 7.659 | 0.0024763712 |
| 736 | 1.884 | -0.059 | 0.019 | -0.108 | -1.493 | 0.213 | NA | 0.925 | 0.338 | NA | 0.971 | 10 | -209.140 | 442.511 | 7.687 | 0.0024411353 |
| 841 | 1.937 | NA | NA | NA | -1.373 | NA | NA | 0.855 | NA | 0.001 | 0.904 | 6 | -214.515 | 442.530 | 7.707 | 0.0024174596 |
| 720 | 1.828 | -0.005 | 0.081 | 0.046 | -1.394 | NA | NA | 0.885 | 0.262 | NA | 1.024 | 9 | -210.597 | 442.589 | 7.766 | 0.0023472717 |
| 604 | 1.898 | -0.011 | 0.213 | NA | -1.363 | 0.127 | NA | 0.845 | NA | NA | 0.930 | 8 | -211.991 | 442.650 | 7.826 | 0.0022775518 |
| 607 | 1.891 | NA | 0.212 | 0.010 | -1.356 | 0.123 | NA | 0.839 | NA | NA | 0.923 | 8 | -211.994 | 442.655 | 7.831 | 0.0022715015 |
| 857 | 1.943 | NA | NA | NA | -1.388 | 0.179 | NA | 0.869 | NA | 0.067 | 0.862 | 7 | -213.324 | 442.684 | 7.860 | 0.0022391394 |
| 988 | 1.843 | -0.022 | -0.005 | NA | -1.467 | 0.157 | NA | 0.931 | 0.351 | -0.073 | 0.995 | 10 | -209.233 | 442.698 | 7.874 | 0.0022234472 |
| 633 | 1.918 | NA | NA | NA | -1.407 | 0.163 | -0.077 | 0.860 | NA | NA | 0.911 | 7 | -213.338 | 442.712 | 7.889 | 0.0022071363 |
| 605 | 1.932 | NA | NA | 0.054 | -1.385 | 0.140 | NA | 0.885 | NA | NA | 0.894 | 7 | -213.370 | 442.777 | 7.953 | 0.0021375320 |
| 602 | 1.953 | -0.024 | NA | NA | -1.406 | 0.162 | NA | 0.895 | NA | NA | 0.896 | 7 | -213.419 | 442.874 | 8.050 | 0.0020356294 |
| 137 | 2.438 | NA | NA | NA | -0.708 | NA | NA | NA | 0.269 | NA | NA | 4 | -217.175 | 443.040 | 8.216 | 0.0018740150 |
| 221 | 2.377 | NA | NA | -0.129 | -1.076 | 0.240 | NA | 0.398 | 0.326 | NA | NA | 7 | -213.558 | 443.152 | 8.328 | 0.0017715904 |
| 875 | 1.837 | NA | 0.261 | NA | -1.348 | NA | -0.148 | 0.782 | NA | 0.093 | 1.027 | 8 | -212.267 | 443.201 | 8.377 | 0.0017288359 |
| 623 | 1.841 | NA | 0.228 | 0.078 | -1.352 | NA | -0.104 | 0.834 | NA | NA | 1.037 | 8 | -212.270 | 443.207 | 8.383 | 0.0017234459 |
| 153 | 2.420 | NA | NA | NA | -0.715 | 0.158 | NA | NA | 0.287 | NA | NA | 5 | -216.122 | 443.298 | 8.474 | 0.0016472706 |
| 218 | 2.307 | 0.077 | NA | NA | -1.021 | 0.193 | NA | 0.376 | 0.292 | NA | NA | 7 | -213.726 | 443.488 | 8.665 | 0.0014974741 |
| 847 | 1.874 | NA | 0.228 | 0.075 | -1.322 | NA | NA | 0.837 | NA | 0.047 | 0.974 | 8 | -212.421 | 443.509 | 8.685 | 0.0014822040 |
| 620 | 1.863 | -0.019 | 0.246 | NA | -1.364 | NA | -0.098 | 0.824 | NA | NA | 1.020 | 8 | -212.452 | 443.570 | 8.747 | 0.0014372034 |
| 592 | 1.863 | 0.019 | 0.223 | 0.082 | -1.320 | NA | NA | 0.833 | NA | NA | 0.960 | 8 | -212.470 | 443.606 | 8.782 | 0.0014118203 |
| 621 | 1.886 | NA | NA | 0.136 | -1.376 | NA | -0.092 | 0.881 | NA | NA | 1.010 | 7 | -213.816 | 443.669 | 8.845 | 0.0013682316 |
| 891 | 1.837 | NA | 0.233 | NA | -1.378 | 0.165 | -0.187 | 0.788 | NA | 0.168 | 1.007 | 9 | -211.214 | 443.823 | 9.000 | 0.0012664062 |
| 844 | 1.891 | -0.012 | 0.245 | NA | -1.333 | NA | NA | 0.824 | NA | 0.038 | 0.957 | 8 | -212.595 | 443.857 | 9.033 | 0.0012453465 |
| 590 | 1.900 | 0.025 | NA | 0.144 | -1.344 | NA | NA | 0.874 | NA | NA | 0.934 | 7 | -213.943 | 443.923 | 9.099 | 0.0012050360 |
| 845 | 1.916 | NA | NA | 0.132 | -1.356 | NA | NA | 0.891 | NA | 0.022 | 0.956 | 7 | -213.952 | 443.940 | 9.116 | 0.0011948853 |
| 473 | 2.395 | NA | NA | NA | -1.021 | 0.176 | NA | 0.395 | 0.304 | -0.041 | NA | 7 | -213.996 | 444.028 | 9.205 | 0.0011430536 |
| 219 | 2.395 | NA | 0.014 | NA | -1.013 | 0.183 | NA | 0.389 | 0.283 | NA | NA | 7 | -214.035 | 444.107 | 9.284 | 0.0010988768 |
| 249 | 2.396 | NA | NA | NA | -1.017 | 0.186 | 0.006 | 0.396 | 0.289 | NA | NA | 7 | -214.039 | 444.115 | 9.291 | 0.0010946917 |
| 457 | 2.416 | NA | NA | NA | -0.995 | NA | NA | 0.361 | 0.312 | -0.119 | NA | 6 | -215.317 | 444.135 | 9.311 | 0.0010839898 |
| 649 | 2.279 | NA | NA | NA | -0.727 | NA | NA | NA | 0.282 | NA | 0.286 | 5 | -216.557 | 444.166 | 9.342 | 0.0010672474 |
| 11 | 2.446 | NA | 0.242 | NA | -0.679 | NA | NA | NA | NA | NA | NA | 4 | -217.863 | 444.416 | 9.592 | 0.0009415836 |
| 203 | 2.417 | NA | 0.094 | NA | -0.948 | NA | NA | 0.321 | 0.224 | NA | NA | 6 | -215.480 | 444.460 | 9.636 | 0.0009211244 |
| 202 | 2.356 | 0.058 | NA | NA | -0.972 | NA | NA | 0.323 | 0.266 | NA | NA | 6 | -215.538 | 444.576 | 9.753 | 0.0008691587 |
| 75 | 2.434 | NA | 0.216 | NA | -0.913 | NA | NA | 0.304 | NA | NA | NA | 5 | -216.763 | 444.578 | 9.754 | 0.0008685717 |
| 139 | 2.430 | NA | 0.131 | NA | -0.698 | NA | NA | NA | 0.213 | NA | NA | 5 | -216.764 | 444.582 | 9.758 | 0.0008668746 |
| 154 | 2.290 | 0.112 | NA | NA | -0.750 | 0.169 | NA | NA | 0.285 | NA | NA | 6 | -215.549 | 444.598 | 9.774 | 0.0008599298 |
| 889 | 1.896 | NA | NA | NA | -1.412 | 0.196 | -0.159 | 0.839 | NA | 0.135 | 0.943 | 8 | -212.972 | 444.610 | 9.787 | 0.0008544811 |
# Model averaging Version 1: use all models with delta AIC score less than 4
model.set.chin.4 <- get.models(model.set.chin, subset = delta < 4)
avg_model.chin.4 <- model.avg(model.set.chin.4)
summary(avg_model.chin.4)
##
## Call:
## model.avg(object = model.set.chin.4)
##
## Component model call:
## glm.nb(formula = abundance ~ <13 unique rhs>, data = b.chin.hab,
## na.action = na.fail, init.theta = <13 unique values>, link = log)
##
## Component models:
## df logLik AICc delta weight
## 4/5/7/8/10 7 -209.39 434.82 0.00 0.21
## 4/7/8/10 6 -210.90 435.31 0.48 0.16
## 4/7/8/9/10 7 -210.39 436.81 1.99 0.08
## 4/5/6/7/8/10 8 -209.08 436.82 2.00 0.08
## 4/5/7/8/9/10 8 -209.26 437.18 2.36 0.06
## 3/4/5/7/8/10 8 -209.28 437.23 2.41 0.06
## 4/6/7/8/10 7 -210.67 437.39 2.56 0.06
## 2/4/7/8/10 7 -210.68 437.40 2.58 0.06
## 2/4/5/7/8/10 8 -209.38 437.42 2.60 0.06
## 1/4/5/7/8/10 8 -209.38 437.43 2.60 0.06
## 3/4/7/8/10 7 -210.76 437.56 2.74 0.05
## 1/4/7/8/10 7 -210.85 437.73 2.90 0.05
## 2/4/7/10 6 -212.64 438.79 3.96 0.03
##
## Term codes:
## j2 meanturb s.do s.J.date s.pH s.sal s.temp stnwidth
## 1 2 3 4 5 6 7 8
## vegelev Year
## 9 10
##
## Model-averaged coefficients:
## (full average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 1.833302 0.201193 0.205535 8.920 < 2e-16 ***
## s.J.date -1.436520 0.249962 0.255335 5.626 < 2e-16 ***
## s.pH 0.090995 0.126463 0.127902 0.711 0.47681
## s.temp 0.895418 0.270309 0.276154 3.242 0.00119 **
## stnwidth 0.310198 0.133734 0.136121 2.279 0.02268 *
## Year 0.990602 0.303815 0.310332 3.192 0.00141 **
## vegelev -0.014297 0.061036 0.061875 0.231 0.81727
## s.sal -0.015460 0.066537 0.067484 0.229 0.81880
## s.do -0.000654 0.048585 0.049450 0.013 0.98945
## meanturb 0.013482 0.067085 0.067905 0.199 0.84262
## j2 -0.002526 0.031215 0.031854 0.079 0.93680
##
## (conditional average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 1.833302 0.201193 0.205535 8.920 < 2e-16 ***
## s.J.date -1.436520 0.249962 0.255335 5.626 < 2e-16 ***
## s.pH 0.175321 0.126608 0.129364 1.355 0.17534
## s.temp 0.895418 0.270309 0.276154 3.242 0.00119 **
## stnwidth 0.319258 0.124558 0.127191 2.510 0.01207 *
## Year 0.990602 0.303815 0.310332 3.192 0.00141 **
## vegelev -0.102398 0.132895 0.135646 0.755 0.45031
## s.sal -0.116291 0.146887 0.150102 0.775 0.43849
## s.do -0.005729 0.143703 0.146264 0.039 0.96876
## meanturb 0.095407 0.155018 0.157524 0.606 0.54474
## j2 -0.024218 0.093901 0.095938 0.252 0.80070
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Relative variable importance:
## s.J.date s.temp Year stnwidth s.pH meanturb vegelev
## Importance: 1.00 1.00 1.00 0.97 0.52 0.14 0.14
## N containing models: 13 13 13 12 6 3 2
## s.sal s.do j2
## Importance: 0.13 0.11 0.10
## N containing models: 2 2 2
chin.ci <- data.frame(confint(avg_model.chin.4, full = TRUE))
# Get pseudo R squared values for models up to delta < 4
chin.4.Rsq <- rsquared(model.set.chin.4, aicc = TRUE)
## write tables to .csv for easy comparison and plugging into appendix table
avg_model_4df.chin <- data.frame(avg_model.chin.4$msTable)
avg_model_components4.chin <- cbind(chin.4.Rsq, avg_model_4df.chin)
r = data.frame(Coeff = rownames(avg_model_4df.chin, rep(NA, length(avg_model_components4.chin))))
avg_model_components4.chin <- cbind(avg_model_components4.chin, r)
avg_model_components4.chin <- avg_model_components4.chin[, -c(6, 7)]
# write.csv(avg_model_components4.chin,
# '/Users/Lia/Documents/Git/Fraser-salmon/all.data/avg_model_components4_beachchin.csv')Parameter Plot
The results of model averaging including all top ranked models up to delta AICc 4
## Warning: package 'cowplot' was built under R version 3.4.3##
## Attaching package: 'cowplot'## The following objects are masked from 'package:sjPlot':
##
## plot_grid, save_plot## The following object is masked from 'package:ggplot2':
##
## ggsave
Check fit of top model by avg – including all parameters with weight >0.5
avg.chin <- glm.nb(abundance ~ s.J.date + Year + s.temp + stnwidth + s.pH, data = b.chin.hab,
na.action = "na.fail")
summary(avg.chin) # AIC 432.9, deviance explained = 100*((129.963-64.196)/129.963) = 50.6 %##
## Call:
## glm.nb(formula = abundance ~ s.J.date + Year + s.temp + stnwidth +
## s.pH, data = b.chin.hab, na.action = "na.fail", init.theta = 1.414591464,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5158 -0.9501 -0.4285 0.3706 1.7588
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.8395 0.1966 9.357 < 2e-16 ***
## s.J.date -1.4488 0.2459 -5.892 0.0000000038 ***
## Year 0.9591 0.2963 3.237 0.001210 **
## s.temp 0.9078 0.2669 3.401 0.000672 ***
## stnwidth 0.3231 0.1185 2.725 0.006427 **
## s.pH 0.1744 0.1245 1.401 0.161140
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.4146) family taken to be 1)
##
## Null deviance: 129.963 on 62 degrees of freedom
## Residual deviance: 64.196 on 57 degrees of freedom
## AIC: 432.79
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.415
## Std. Err.: 0.270
##
## 2 x log-likelihood: -418.787vif(avg.chin) # - Jday and temp still above 5 (5.055 and 5.93, respectively)## s.J.date Year s.temp stnwidth s.pH
## 5.055297 1.639788 5.930079 1.067749 1.050073### Plot residuals vs. fitted values
plot(fitted(avg.chin), resid(avg.chin), main = "Averaged Chinook GLM", xlab = "Fitted",
ylab = "Pearson residuals")
## q plot
qqnorm(x = b.chin.hab$abundance, y = resid(avg.chin), main = "Q-Q Averaged Chinook GLM")
qqline(resid(avg.chin), col = 2)
Chum beach seine model
response = average abundance of chum salmon [per beach seine site]
This is 2nd of 4 beach models for (1)Chinook, (2)chum, (3)other migratory species and (4)resident species
Load all data - standardize variables, create subgroups from beach seine data
### Separate into groups; automatically removes 3 x unidentified larval and 28
### x '0' abundance rows 2
b.2 <- b[which(b$Species == "Chum"), ]
### standardize continuous variables
b.2$s.temp <- standardize(b.2$Temp.surf, centerFun = mean, scaleFun = sd)
b.2$s.sal <- standardize(b.2$Sal.surf, centerFun = mean, scaleFun = sd)
b.2$s.do <- standardize(b.2$DOmg.surf, centerFun = mean, scaleFun = sd)
b.2$s.pH <- standardize(b.2$pH.surf, centerFun = mean, scaleFun = sd)
b.2$s.J.date <- standardize(b.2$J.date, centerFun = mean, scaleFun = sd)
### Create variable j2 which is Julian day squared
b.2$j2 <- b.2$s.J.date^2
# 1 summarize by site-day
b.chum <- ddply(b.2, .(Year, J.date, s.J.date, j2, Habitat, Site, s.temp, s.sal,
s.do, s.pH), summarize, abundance = sum(abundance))
# plot abundance over julian day to see distribution:
plot(abundance ~ J.date, data = b.chum)
Habitat variables
## Add marsh habitat variables for each site
b.chum.hab <- b.chum
b.chum.hab$stnwidth <- hab2[match(b.chum.hab$Site, hab$Site), 2]
b.chum.hab$vegelev <- hab2[match(b.chum.hab$Site, hab$Site), 3]
b.chum.hab$shtdensity <- hab2[match(b.chum.hab$Site, hab$Site), 4]
b.chum.hab$shtdenhigh <- hab2[match(b.chum.hab$Site, hab$Site), 5]
b.chum.hab$tcelev <- hab2[match(b.chum.hab$Site, hab$Site), 6]
b.chum.hab$angbank <- hab2[match(b.chum.hab$Site, hab$Site), 7]
b.chum.hab$meanturb <- hab2[match(b.chum.hab$Site, hab$Site), 8]
## standardize habitat variables using 'standardize' function from robustHD
## package
b.chum.hab$stnwidth <- standardize(b.chum.hab$stnwidth, centerFun = mean, scaleFun = sd)
b.chum.hab$stnwidth <- as.numeric(b.chum.hab$stnwidth)
b.chum.hab$vegelev <- standardize(b.chum.hab$vegelev, centerFun = mean, scaleFun = sd)
b.chum.hab$vegelev <- as.numeric(b.chum.hab$vegelev)
b.chum.hab$shtdensity <- standardize(b.chum.hab$shtdensity, centerFun = mean,
scaleFun = sd)
b.chum.hab$shtdensity <- as.numeric(b.chum.hab$shtdensity)
b.chum.hab$shtdenhigh <- standardize(b.chum.hab$shtdenhigh, centerFun = mean,
scaleFun = sd)
b.chum.hab$shtdenhigh <- as.numeric(b.chum.hab$shtdenhigh)
b.chum.hab$tcelev <- standardize(b.chum.hab$tcelev, centerFun = mean, scaleFun = sd)
b.chum.hab$tcelev <- as.numeric(b.chum.hab$tcelev)
b.chum.hab$angbank <- standardize(b.chum.hab$angbank, centerFun = mean, scaleFun = sd)
b.chum.hab$angbank <- as.numeric(b.chum.hab$angbank)
b.chum.hab$meanturb <- standardize(b.chum.hab$meanturb, centerFun = mean, scaleFun = sd)
b.chum.hab$meanturb <- as.numeric(b.chum.hab$meanturb)VIF for collinearity of habitat variables
Assess variance inflation factors
# alias function identifies covariates that are multiples of each other - in
# this case have some habitat parameters causing issues.
alias(lm(abundance ~ s.J.date + j2 + Year + s.temp + s.sal + s.do + s.pH + stnwidth +
vegelev + shtdenhigh + shtdensity + tcelev + angbank + meanturb, data = b.chum.hab))
## Model :
## abundance ~ s.J.date + j2 + Year + s.temp + s.sal + s.do + s.pH +
## stnwidth + vegelev + shtdenhigh + shtdensity + tcelev + angbank +
## meanturb
##
## Complete :
## (Intercept) s.J.date j2
## tcelev 0 0 0
## angbank 0 0 0
## meanturb 0 0 0
## Year s.temp s.sal
## tcelev 0 0 0
## angbank 0 0 0
## meanturb 0 0 0
## s.do s.pH stnwidth
## tcelev 0 0 8871/14840
## angbank 0 0 -133424/254871
## meanturb 0 0 27432/68227
## vegelev shtdenhigh shtdensity
## tcelev 605155061/1293820358 -376609/756041 1235/3802
## angbank 2957/4166 -9279/22789 -3094/47891
## meanturb -3080706/3772433 -871961/1287673 -35732/62255
# We determined hab variables for Chum should be reduced to meanturb, and
# vegelev based on significance and collinearity when doing model
# exploration.
vif(lm(abundance ~ s.J.date + j2 + Year + s.temp + s.sal + s.do + s.pH + vegelev +
meanturb, data = b.chum.hab)) ##all less than 5
## s.J.date j2 Year s.temp s.sal s.do s.pH vegelev
## 2.603210 2.909987 2.583522 3.466341 3.020429 3.369835 1.630916 2.041842
## meanturb
## 1.165110
## Pearson Corr with all vars
year <- b.chum.hab$Year
jday <- b.chum.hab$s.J.date
j2 <- b.chum.hab$j2
temp <- b.chum.hab$s.temp
do <- b.chum.hab$s.do
ph <- b.chum.hab$s.pH
sal <- b.chum.hab$s.sal
veg <- b.chum.hab$vegelev
turb <- b.chum.hab$meanturb
width <- b.chum.hab$stnwidth
shootden <- b.chum.hab$shtdensity
shtdenhigh <- b.chum.hab$shtdenhigh
elev <- b.chum.hab$tcelev
angbank <- b.chum.hab$angbank
habcovar <- cbind.data.frame(year, jday, j2, temp, do, ph, sal, veg, turb, width,
shootden, shtdenhigh, elev, angbank)
ggpairs(data = na.omit(habcovar), title = "Pearson Correlation plot for all variables in beach seine chum salmon dataset")
ggsave("/Users/Lia/Documents/Git/Fraser-salmon/Ch1.Seasonal_diversity/expl.plots/PearsonCorrAllVars_chummarsh.pdf")Model selection
Full model
## Chum: full model with top habitat variables
chum <- glm.nb(abundance ~ s.J.date + j2 + Year + s.temp + s.sal + s.do + s.pH +
vegelev + meanturb, data = b.chum.hab, na.action = "na.fail")
summary(chum) #AIC 233.64##
## Call:
## glm.nb(formula = abundance ~ s.J.date + j2 + Year + s.temp +
## s.sal + s.do + s.pH + vegelev + meanturb, data = b.chum.hab,
## na.action = "na.fail", init.theta = 0.9805681787, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8251 -1.2597 -0.3236 0.5022 1.4294
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.2029 0.4342 5.073 0.000000391 ***
## s.J.date -0.5285 0.2968 -1.781 0.0749 .
## j2 -0.1948 0.2205 -0.884 0.3770
## Year 0.9408 0.6445 1.460 0.1443
## s.temp 0.1084 0.3458 0.314 0.7539
## s.sal -0.4729 0.3153 -1.500 0.1336
## s.do 0.5665 0.3698 1.532 0.1256
## s.pH 0.3145 0.2707 1.162 0.2452
## vegelev 0.3977 0.2829 1.406 0.1598
## meanturb -0.3448 0.2157 -1.599 0.1098
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.9806) family taken to be 1)
##
## Null deviance: 58.891 on 30 degrees of freedom
## Residual deviance: 33.073 on 21 degrees of freedom
## AIC: 233.64
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.981
## Std. Err.: 0.255
##
## 2 x log-likelihood: -211.642## Now look at the variance inflation factors of this full model
vif(chum)## s.J.date j2 Year s.temp s.sal s.do s.pH vegelev
## 2.366396 2.638545 2.607406 3.100914 3.139654 3.062227 1.634527 2.029153
## meanturb
## 1.210263# Check model assumptions - note that we are not actually assuming linear
# relationships with Jday, but can't tell this general function that.
sjp.glm(chum, type = "ma")## `sjp.glm()` will become deprecated in the future. Please use `plot_model()` instead.## No outliers detected.## NULL
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 0.12753## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.92707## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 2.3502e-16## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at 0.12753## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.92707## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 2.3502e-16## --------------------
## Check significance of terms when they entered the model...## Anova:## Warning in anova.negbin(logreg, test = "Chisq"): tests made without re-
## estimating 'theta'
## Analysis of Deviance Table
##
## Model: Negative Binomial(0.9806), link: log
##
## Response: abundance
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 30 58.891
## s.J.date 1 0.9566 29 57.935 0.3280324
## j2 1 14.0090 28 43.926 0.0001819 ***
## Year 1 3.6769 27 40.249 0.0551719 .
## s.temp 1 0.0207 26 40.228 0.8855816
## s.sal 1 0.5239 25 39.704 0.4691742
## s.do 1 2.8477 24 36.856 0.0915062 .
## s.pH 1 0.2736 23 36.583 0.6009037
## vegelev 1 1.0859 22 35.497 0.2973723
## meanturb 1 2.4243 21 33.073 0.1194640
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# See pseudo R2 for full model
rsquared(chum, aicc = TRUE)## Class Family Link n R.squared
## 1 negbin Negative Binomial log 31 0.4384135Model selection - dredge
Then continued to dredge function to determine optimal model using AICc (chum has 31 observations). Note that due to the small sample size, even the top model has a very low weight (0.08), and consequently many models (42) are included in the averaging set due to very small changes in delta AICc with different iterations.
# Generate model set
model.set.chum <- dredge(chum)
# Create model selection table
model_table.chum <- model.sel(model.set.chum)
options(scipen = 7)
names(model_table.chum) <- c("(Int)", "Jday^2", "Mean turb", "DO", "Jday", "pH",
"Sal", "Temp", "Veg. elev", "Year", "df", "LL", "AICc", "delta", "weight")
kable(head(model_table.chum, n = 100), digits = 3)| (Int) | Jday^2 | Mean turb | DO | Jday | pH | Sal | Temp | Veg. elev | Year | df | LL | AICc | delta | weight | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 13 | 2.621 | NA | NA | 0.923 | -0.367 | NA | NA | NA | NA | NA | 4 | -108.612 | 226.763 | 0.000 | 0.046856556 |
| 5 | 2.646 | NA | NA | 0.838 | NA | NA | NA | NA | NA | NA | 3 | -110.013 | 226.914 | 0.151 | 0.043447042 |
| 69 | 2.599 | NA | NA | 0.886 | NA | NA | NA | -0.319 | NA | NA | 4 | -108.998 | 227.534 | 0.771 | 0.031868966 |
| 15 | 2.572 | NA | -0.302 | 0.917 | -0.377 | NA | NA | NA | NA | NA | 5 | -107.587 | 227.574 | 0.811 | 0.031233524 |
| 7 | 2.610 | NA | -0.290 | 0.847 | NA | NA | NA | NA | NA | NA | 4 | -109.159 | 227.857 | 1.094 | 0.027115205 |
| 269 | 2.306 | NA | NA | 0.900 | -0.411 | NA | NA | NA | NA | 0.474 | 5 | -108.047 | 228.495 | 1.732 | 0.019709864 |
| 133 | 2.635 | NA | NA | 0.885 | NA | NA | NA | NA | 0.201 | NA | 4 | -109.623 | 228.784 | 2.021 | 0.017059660 |
| 37 | 2.610 | NA | NA | 0.831 | NA | NA | 0.189 | NA | NA | NA | 4 | -109.686 | 228.911 | 2.149 | 0.016003777 |
| 77 | 2.592 | NA | NA | 0.920 | -0.280 | NA | NA | -0.207 | NA | NA | 5 | -108.298 | 228.996 | 2.233 | 0.015341997 |
| 261 | 2.433 | NA | NA | 0.826 | NA | NA | NA | NA | NA | 0.329 | 4 | -109.739 | 229.017 | 2.254 | 0.015183630 |
| 6 | 2.806 | -0.173 | NA | 0.629 | NA | NA | NA | NA | NA | NA | 4 | -109.785 | 229.108 | 2.345 | 0.014503333 |
| 2 | 3.170 | -0.516 | NA | NA | NA | NA | NA | NA | NA | NA | 3 | -111.114 | 229.117 | 2.354 | 0.014440692 |
| 71 | 2.575 | NA | -0.240 | 0.881 | NA | NA | NA | -0.269 | NA | NA | 5 | -108.423 | 229.245 | 2.482 | 0.013543770 |
| 141 | 2.611 | NA | NA | 0.944 | -0.340 | NA | NA | NA | 0.139 | NA | 5 | -108.423 | 229.246 | 2.484 | 0.013535084 |
| 18 | 3.109 | -0.473 | NA | NA | NA | 0.313 | NA | NA | NA | NA | 4 | -109.861 | 229.261 | 2.498 | 0.013439000 |
| 66 | 3.116 | -0.540 | NA | NA | NA | NA | NA | -0.376 | NA | NA | 4 | -109.873 | 229.283 | 2.521 | 0.013286939 |
| 14 | 2.735 | -0.128 | NA | 0.767 | -0.352 | NA | NA | NA | NA | NA | 5 | -108.466 | 229.332 | 2.569 | 0.012967934 |
| 21 | 2.649 | NA | NA | 0.796 | NA | 0.075 | NA | NA | NA | NA | 4 | -109.958 | 229.455 | 2.692 | 0.012196116 |
| 45 | 2.603 | NA | NA | 0.908 | -0.342 | NA | 0.091 | NA | NA | NA | 5 | -108.540 | 229.480 | 2.717 | 0.012043351 |
| 70 | 2.796 | -0.225 | NA | 0.612 | NA | NA | NA | -0.345 | NA | NA | 5 | -108.576 | 229.551 | 2.788 | 0.011623464 |
| 29 | 2.623 | NA | NA | 0.894 | -0.365 | 0.046 | NA | NA | NA | NA | 5 | -108.588 | 229.576 | 2.814 | 0.011476432 |
| 135 | 2.595 | NA | -0.303 | 0.901 | NA | NA | NA | NA | 0.221 | NA | 5 | -108.647 | 229.694 | 2.931 | 0.010821136 |
| 197 | 2.589 | NA | NA | 0.923 | NA | NA | NA | -0.308 | 0.173 | NA | 5 | -108.672 | 229.743 | 2.980 | 0.010558866 |
| 271 | 2.301 | NA | -0.289 | 0.896 | -0.417 | NA | NA | NA | NA | 0.416 | 6 | -107.127 | 229.754 | 2.991 | 0.010504088 |
| 258 | 2.733 | -0.522 | NA | NA | NA | NA | NA | NA | NA | 0.654 | 4 | -110.114 | 229.767 | 3.004 | 0.010435365 |
| 20 | 3.045 | -0.446 | -0.340 | NA | NA | 0.440 | NA | NA | NA | NA | 5 | -108.724 | 229.847 | 3.084 | 0.010024106 |
| 82 | 3.070 | -0.501 | NA | NA | NA | 0.292 | NA | -0.336 | NA | NA | 5 | -108.759 | 229.917 | 3.154 | 0.009678037 |
| 39 | 2.572 | NA | -0.292 | 0.842 | NA | NA | 0.193 | NA | NA | NA | 5 | -108.785 | 229.971 | 3.208 | 0.009423040 |
| 143 | 2.559 | NA | -0.306 | 0.945 | -0.347 | NA | NA | NA | 0.153 | NA | 6 | -107.340 | 230.180 | 3.417 | 0.008486253 |
| 23 | 2.611 | NA | -0.332 | 0.746 | NA | 0.186 | NA | NA | NA | NA | 5 | -108.898 | 230.196 | 3.433 | 0.008419450 |
| 31 | 2.574 | NA | -0.340 | 0.821 | -0.374 | 0.157 | NA | NA | NA | NA | 6 | -107.367 | 230.234 | 3.471 | 0.008260439 |
| 325 | 2.685 | NA | NA | 0.898 | NA | NA | NA | -0.368 | NA | -0.144 | 5 | -108.967 | 230.334 | 3.572 | 0.007856236 |
| 263 | 2.440 | NA | -0.279 | 0.838 | NA | NA | NA | NA | NA | 0.264 | 5 | -108.975 | 230.350 | 3.587 | 0.007796556 |
| 85 | 2.601 | NA | NA | 0.868 | NA | 0.031 | NA | -0.315 | NA | NA | 5 | -108.989 | 230.378 | 3.615 | 0.007687679 |
| 101 | 2.596 | NA | NA | 0.883 | NA | NA | 0.027 | -0.305 | NA | NA | 5 | -108.993 | 230.385 | 3.623 | 0.007658219 |
| 79 | 2.558 | NA | -0.281 | 0.914 | -0.323 | NA | NA | -0.133 | NA | NA | 6 | -107.453 | 230.405 | 3.642 | 0.007582632 |
| 8 | 2.734 | -0.133 | -0.276 | 0.688 | NA | NA | NA | NA | NA | NA | 5 | -109.024 | 230.447 | 3.684 | 0.007425643 |
| 47 | 2.553 | NA | -0.302 | 0.904 | -0.352 | NA | 0.095 | NA | NA | NA | 6 | -107.499 | 230.498 | 3.735 | 0.007239276 |
| 266 | 2.603 | -0.523 | NA | NA | -0.326 | NA | NA | NA | NA | 0.802 | 5 | -109.087 | 230.573 | 3.811 | 0.006971391 |
| 26 | 3.097 | -0.485 | NA | NA | -0.269 | 0.316 | NA | NA | NA | NA | 5 | -109.096 | 230.591 | 3.828 | 0.006909639 |
| 16 | 2.634 | -0.069 | -0.291 | 0.836 | -0.367 | NA | NA | NA | NA | NA | 6 | -107.546 | 230.593 | 3.830 | 0.006903870 |
| 10 | 3.161 | -0.526 | NA | NA | -0.242 | NA | NA | NA | NA | NA | 4 | -110.546 | 230.631 | 3.868 | 0.006773211 |
| 4 | 3.142 | -0.512 | -0.219 | NA | NA | NA | NA | NA | NA | NA | 4 | -110.648 | 230.834 | 4.071 | 0.006119945 |
| 34 | 3.101 | -0.498 | NA | NA | NA | NA | 0.236 | NA | NA | NA | 4 | -110.711 | 230.960 | 4.197 | 0.005745496 |
| 270 | 2.423 | -0.194 | NA | 0.658 | -0.389 | NA | NA | NA | NA | 0.554 | 6 | -107.731 | 230.961 | 4.199 | 0.005742065 |
| 262 | 2.581 | -0.244 | NA | 0.522 | NA | NA | NA | NA | NA | 0.446 | 5 | -109.324 | 231.047 | 4.284 | 0.005500541 |
| 274 | 2.808 | -0.484 | NA | NA | NA | 0.265 | NA | NA | NA | 0.467 | 5 | -109.340 | 231.081 | 4.318 | 0.005409451 |
| 389 | 2.454 | NA | NA | 0.870 | NA | NA | NA | NA | 0.178 | 0.281 | 5 | -109.424 | 231.249 | 4.486 | 0.004973398 |
| 28 | 3.015 | -0.450 | -0.342 | NA | -0.280 | 0.438 | NA | NA | NA | NA | 6 | -107.875 | 231.251 | 4.488 | 0.004968677 |
| 38 | 2.770 | -0.177 | NA | 0.615 | NA | NA | 0.195 | NA | NA | NA | 5 | -109.438 | 231.277 | 4.514 | 0.004904452 |
| 149 | 2.637 | NA | NA | 0.817 | NA | 0.124 | NA | NA | 0.230 | NA | 5 | -109.473 | 231.347 | 4.584 | 0.004735403 |
| 397 | 2.318 | NA | NA | 0.919 | -0.389 | NA | NA | NA | 0.105 | 0.445 | 6 | -107.934 | 231.367 | 4.605 | 0.004686947 |
| 22 | 2.903 | -0.272 | NA | 0.404 | NA | 0.186 | NA | NA | NA | NA | 5 | -109.514 | 231.428 | 4.666 | 0.004546139 |
| 134 | 2.736 | -0.108 | NA | 0.747 | NA | NA | NA | NA | 0.169 | NA | 5 | -109.539 | 231.478 | 4.715 | 0.004435294 |
| 301 | 2.247 | NA | NA | 0.909 | -0.447 | NA | -0.099 | NA | NA | 0.591 | 6 | -107.992 | 231.483 | 4.720 | 0.004423265 |
| 199 | 2.563 | NA | -0.249 | 0.922 | NA | NA | NA | -0.248 | 0.187 | NA | 6 | -108.027 | 231.554 | 4.791 | 0.004270692 |
| 165 | 2.620 | NA | NA | 0.867 | NA | NA | 0.093 | NA | 0.140 | NA | 5 | -109.579 | 231.558 | 4.795 | 0.004261933 |
| 333 | 2.278 | NA | NA | 0.899 | -0.431 | NA | NA | 0.038 | NA | 0.524 | 6 | -108.043 | 231.587 | 4.824 | 0.004199797 |
| 285 | 2.309 | NA | NA | 0.895 | -0.410 | 0.008 | NA | NA | NA | 0.471 | 6 | -108.047 | 231.594 | 4.831 | 0.004185970 |
| 78 | 2.732 | -0.163 | NA | 0.721 | -0.243 | NA | NA | -0.234 | NA | NA | 6 | -108.071 | 231.643 | 4.880 | 0.004084400 |
| 293 | 2.502 | NA | NA | 0.828 | NA | NA | 0.134 | NA | NA | 0.183 | 5 | -109.626 | 231.652 | 4.889 | 0.004065417 |
| 68 | 3.098 | -0.533 | -0.150 | NA | NA | NA | NA | -0.340 | NA | NA | 5 | -109.658 | 231.716 | 4.953 | 0.003937852 |
| 322 | 2.894 | -0.538 | NA | NA | NA | NA | NA | -0.265 | NA | 0.361 | 5 | -109.663 | 231.726 | 4.963 | 0.003917347 |
| 84 | 3.026 | -0.472 | -0.261 | NA | NA | 0.382 | NA | -0.255 | NA | NA | 6 | -108.117 | 231.734 | 4.971 | 0.003902126 |
| 50 | 3.071 | -0.465 | NA | NA | NA | 0.292 | 0.152 | NA | NA | NA | 5 | -109.669 | 231.737 | 4.974 | 0.003895890 |
| 53 | 2.613 | NA | NA | 0.807 | NA | 0.043 | 0.180 | NA | NA | NA | 5 | -109.669 | 231.738 | 4.975 | 0.003894211 |
| 205 | 2.583 | NA | NA | 0.942 | -0.253 | NA | NA | -0.204 | 0.132 | NA | 6 | -108.120 | 231.739 | 4.976 | 0.003891784 |
| 130 | 3.170 | -0.516 | NA | NA | NA | NA | NA | NA | -0.022 | NA | 4 | -111.109 | 231.757 | 4.994 | 0.003856916 |
| 72 | 2.740 | -0.184 | -0.210 | 0.659 | NA | NA | NA | -0.294 | NA | NA | 6 | -108.148 | 231.796 | 5.033 | 0.003783293 |
| 151 | 2.591 | NA | -0.353 | 0.767 | NA | 0.246 | NA | NA | 0.265 | NA | 6 | -108.170 | 231.840 | 5.077 | 0.003700746 |
| 277 | 2.446 | NA | NA | 0.802 | NA | 0.043 | NA | NA | NA | 0.310 | 5 | -109.723 | 231.846 | 5.083 | 0.003690078 |
| 146 | 3.094 | -0.462 | NA | NA | NA | 0.345 | NA | NA | 0.127 | NA | 5 | -109.724 | 231.848 | 5.085 | 0.003686308 |
| 260 | 2.735 | -0.520 | -0.199 | NA | NA | NA | NA | NA | NA | 0.616 | 5 | -109.736 | 231.873 | 5.110 | 0.003640283 |
| 282 | 2.677 | -0.489 | NA | NA | -0.336 | 0.263 | NA | NA | NA | 0.620 | 6 | -108.192 | 231.883 | 5.120 | 0.003621802 |
| 74 | 3.111 | -0.537 | NA | NA | -0.091 | NA | NA | -0.337 | NA | NA | 5 | -109.806 | 232.013 | 5.250 | 0.003394477 |
| 93 | 2.594 | NA | NA | 0.898 | -0.281 | 0.035 | NA | -0.203 | NA | NA | 6 | -108.285 | 232.070 | 5.307 | 0.003299132 |
| 194 | 3.120 | -0.543 | NA | NA | NA | NA | NA | -0.378 | -0.049 | NA | 5 | -109.847 | 232.095 | 5.332 | 0.003257983 |
| 109 | 2.592 | NA | NA | 0.919 | -0.280 | NA | 0.004 | -0.205 | NA | NA | 6 | -108.298 | 232.096 | 5.333 | 0.003256704 |
| 30 | 2.803 | -0.199 | NA | 0.600 | -0.337 | 0.131 | NA | NA | NA | NA | 6 | -108.316 | 232.133 | 5.370 | 0.003196857 |
| 87 | 2.579 | NA | -0.268 | 0.811 | NA | 0.118 | NA | -0.245 | NA | NA | 6 | -108.318 | 232.137 | 5.374 | 0.003190311 |
| 98 | 3.112 | -0.538 | NA | NA | NA | NA | 0.018 | -0.367 | NA | NA | 5 | -109.871 | 232.142 | 5.379 | 0.003182667 |
| 157 | 2.613 | NA | NA | 0.889 | -0.333 | 0.093 | NA | NA | 0.170 | NA | 6 | -108.332 | 232.165 | 5.402 | 0.003146370 |
| 142 | 2.691 | -0.089 | NA | 0.832 | -0.335 | NA | NA | NA | 0.111 | NA | 6 | -108.360 | 232.219 | 5.456 | 0.003061326 |
| 86 | 2.873 | -0.303 | NA | 0.420 | NA | 0.159 | NA | -0.333 | NA | NA | 6 | -108.380 | 232.260 | 5.497 | 0.003000189 |
| 327 | 2.666 | NA | -0.239 | 0.892 | NA | NA | NA | -0.321 | NA | -0.152 | 6 | -108.388 | 232.277 | 5.514 | 0.002975030 |
| 46 | 2.717 | -0.129 | NA | 0.750 | -0.326 | NA | 0.093 | NA | NA | NA | 6 | -108.391 | 232.283 | 5.520 | 0.002966136 |
| 103 | 2.567 | NA | -0.246 | 0.875 | NA | NA | 0.065 | -0.233 | NA | NA | 6 | -108.393 | 232.286 | 5.523 | 0.002961601 |
| 173 | 2.612 | NA | NA | 0.946 | -0.341 | NA | -0.004 | NA | 0.142 | NA | 6 | -108.423 | 232.346 | 5.583 | 0.002873089 |
| 148 | 3.017 | -0.429 | -0.356 | NA | NA | 0.489 | NA | NA | 0.172 | NA | 6 | -108.437 | 232.373 | 5.610 | 0.002834510 |
| 198 | 2.744 | -0.176 | NA | 0.697 | NA | NA | NA | -0.332 | 0.119 | NA | 6 | -108.438 | 232.377 | 5.614 | 0.002829707 |
| 229 | 2.612 | NA | NA | 0.988 | NA | NA | -0.259 | -0.425 | 0.318 | NA | 6 | -108.439 | 232.378 | 5.615 | 0.002828478 |
| 24 | 2.864 | -0.269 | -0.330 | 0.359 | NA | 0.307 | NA | NA | NA | NA | 6 | -108.441 | 232.381 | 5.619 | 0.002822976 |
| 276 | 2.845 | -0.458 | -0.304 | NA | NA | 0.383 | NA | NA | NA | 0.323 | 6 | -108.473 | 232.446 | 5.683 | 0.002733817 |
| 61 | 2.606 | NA | NA | 0.885 | -0.343 | 0.038 | 0.086 | NA | NA | NA | 6 | -108.524 | 232.548 | 5.785 | 0.002597375 |
| 391 | 2.461 | NA | -0.294 | 0.891 | NA | NA | NA | NA | 0.205 | 0.211 | 6 | -108.529 | 232.558 | 5.795 | 0.002584493 |
| 12 | 3.123 | -0.518 | -0.214 | NA | -0.242 | NA | NA | NA | NA | NA | 5 | -110.083 | 232.567 | 5.804 | 0.002573197 |
| 386 | 2.733 | -0.523 | NA | NA | NA | NA | NA | NA | -0.031 | 0.655 | 5 | -110.104 | 232.608 | 5.845 | 0.002520464 |
| 90 | 3.060 | -0.495 | NA | NA | -0.152 | 0.301 | NA | -0.272 | NA | NA | 6 | -108.557 | 232.615 | 5.852 | 0.002512494 |
| 290 | 2.745 | -0.520 | NA | NA | NA | NA | 0.026 | NA | NA | 0.625 | 5 | -110.111 | 232.622 | 5.859 | 0.002503511 |
| 102 | 2.795 | -0.225 | NA | 0.612 | NA | NA | 0.006 | -0.342 | NA | NA | 6 | -108.575 | 232.651 | 5.888 | 0.002467574 |
# Model averaging Version 1: use all models with delta AIC score less than 4
model.set.chum.4 <- get.models(model.set.chum, subset = delta < 4)
avg_model.chum.4 <- model.avg(model.set.chum.4)
summary(avg_model.chum.4)
##
## Call:
## model.avg(object = model.set.chum.4)
##
## Component model call:
## glm.nb(formula = abundance ~ <42 unique rhs>, data = b.chum.hab,
## na.action = na.fail, init.theta = <42 unique values>, link = log)
##
## Component models:
## df logLik AICc delta weight
## 34 4 -108.61 226.76 0.00 0.08
## 3 3 -110.01 226.91 0.15 0.07
## 37 4 -109.00 227.53 0.77 0.05
## 234 5 -107.59 227.57 0.81 0.05
## 23 4 -109.16 227.86 1.09 0.05
## 349 5 -108.05 228.49 1.73 0.03
## 38 4 -109.62 228.78 2.02 0.03
## 36 4 -109.69 228.91 2.15 0.03
## 347 5 -108.30 229.00 2.23 0.03
## 39 4 -109.74 229.02 2.25 0.03
## 13 4 -109.78 229.11 2.35 0.02
## 1 3 -111.11 229.12 2.35 0.02
## 237 5 -108.42 229.25 2.48 0.02
## 348 5 -108.42 229.25 2.48 0.02
## 15 4 -109.86 229.26 2.50 0.02
## 17 4 -109.87 229.28 2.52 0.02
## 134 5 -108.47 229.33 2.57 0.02
## 35 4 -109.96 229.45 2.69 0.02
## 346 5 -108.54 229.48 2.72 0.02
## 137 5 -108.58 229.55 2.79 0.02
## 345 5 -108.59 229.58 2.81 0.02
## 238 5 -108.65 229.69 2.93 0.02
## 378 5 -108.67 229.74 2.98 0.02
## 2349 6 -107.13 229.75 2.99 0.02
## 19 4 -110.11 229.77 3.00 0.02
## 125 5 -108.72 229.85 3.08 0.02
## 157 5 -108.76 229.92 3.15 0.02
## 236 5 -108.79 229.97 3.21 0.02
## 2348 6 -107.34 230.18 3.42 0.01
## 235 5 -108.90 230.20 3.43 0.01
## 2345 6 -107.37 230.23 3.47 0.01
## 379 5 -108.97 230.33 3.57 0.01
## 239 5 -108.97 230.35 3.59 0.01
## 357 5 -108.99 230.38 3.61 0.01
## 367 5 -108.99 230.39 3.62 0.01
## 2347 6 -107.45 230.41 3.64 0.01
## 123 5 -109.02 230.45 3.68 0.01
## 2346 6 -107.50 230.50 3.74 0.01
## 149 5 -109.09 230.57 3.81 0.01
## 145 5 -109.10 230.59 3.83 0.01
## 1234 6 -107.55 230.59 3.83 0.01
## 14 4 -110.55 230.63 3.87 0.01
##
## Term codes:
## j2 meanturb s.do s.J.date s.pH s.sal s.temp vegelev
## 1 2 3 4 5 6 7 8
## Year
## 9
##
## Model-averaged coefficients:
## (full average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 2.66487 0.31503 0.32401 8.225 <2e-16 ***
## s.do 0.72465 0.39016 0.39608 1.830 0.0673 .
## s.J.date -0.14133 0.22008 0.22390 0.631 0.5279
## s.temp -0.07041 0.17154 0.17492 0.403 0.6873
## meanturb -0.08859 0.17905 0.18265 0.485 0.6277
## Year 0.05483 0.23053 0.23587 0.232 0.8162
## vegelev 0.01855 0.08886 0.09141 0.203 0.8392
## s.sal 0.01168 0.07314 0.07546 0.155 0.8770
## j2 -0.09338 0.20253 0.20422 0.457 0.6475
## s.pH 0.03107 0.13235 0.13566 0.229 0.8188
##
## (conditional average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 2.6649 0.3150 0.3240 8.225 < 2e-16 ***
## s.do 0.8600 0.2535 0.2642 3.255 0.00113 **
## s.J.date -0.3547 0.2142 0.2239 1.584 0.11314
## s.temp -0.3010 0.2375 0.2478 1.215 0.22445
## meanturb -0.2961 0.2138 0.2237 1.324 0.18561
## Year 0.4084 0.5014 0.5196 0.786 0.43190
## vegelev 0.1792 0.2179 0.2279 0.786 0.43162
## s.sal 0.1303 0.2103 0.2193 0.594 0.55234
## j2 -0.3753 0.2430 0.2486 1.509 0.13117
## s.pH 0.2061 0.2831 0.2933 0.703 0.48223
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Relative variable importance:
## s.do s.J.date meanturb j2 s.temp s.pH Year vegelev
## Importance: 0.84 0.40 0.30 0.25 0.23 0.15 0.13 0.10
## N containing models: 33 17 15 14 11 9 7 5
## s.sal
## Importance: 0.09
## N containing models: 5
chum.ci <- data.frame(confint(avg_model.chum.4, full = TRUE))
# Get pseudo R squared values for models up to delta < 4
chum.4.Rsq <- rsquared(model.set.chum.4, aicc = TRUE)
## write tables to .csv for easy comparison and plugging into appendix table
avg_model_4df.chum <- data.frame(avg_model.chum.4$msTable)
avg_model_components4.chum <- cbind(chum.4.Rsq, avg_model_4df.chum)
r = data.frame(Coeff = rownames(avg_model_4df.chum, rep(NA, length(avg_model_components4.chum))))
avg_model_components4.chum <- cbind(avg_model_components4.chum, r)
avg_model_components4.chum <- avg_model_components4.chum[, -c(6, 7)]
# write.csv(avg_model_components4.chum,
# '/Users/Lia/Documents/Git/Fraser-salmon/all.data/avg_model_components4_beachchum.csv')Parameter Plot
The results of model averaging including all top ranked models up to delta AICc 4
Check fit of top model by avg – including all parameters with weight >0.5
avg.chum <- glm.nb(abundance ~ s.do, data = b.chum.hab, na.action = "na.fail")
summary(avg.chum) # AIC 226.03 ##
## Call:
## glm.nb(formula = abundance ~ s.do, data = b.chum.hab, na.action = "na.fail",
## init.theta = 0.7611315465, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7338 -1.2327 -0.4714 0.4236 1.7209
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.6465 0.2152 12.298 < 2e-16 ***
## s.do 0.8380 0.2296 3.651 0.000262 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.7611) family taken to be 1)
##
## Null deviance: 46.900 on 30 degrees of freedom
## Residual deviance: 34.031 on 29 degrees of freedom
## AIC: 226.03
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.761
## Std. Err.: 0.186
##
## 2 x log-likelihood: -220.025# vif(avg.chum) ## Can't run VIF because only one variable!
### Plot residuals vs. fitted values
plot(fitted(avg.chum), resid(avg.chum), main = "Averaged Chum GLM", xlab = "Fitted",
ylab = "Pearson residuals")
## q plot
qqnorm(x = b.chum.hab$abundance, y = resid(avg.chum), main = "Q-Q Averaged Chum GLM")
qqline(resid(avg.chum), col = 2)
Migratory beach seine model
response = average abundance of migratory fishes [per beach seine site]
This is 3rd of 4 beach models for (1)Chinook, (2)migra, (3)other migratory species and (4)resident species
Load all data - standardize variables, create subgroups from beach seine data
### Grab just migratory catch - automatically removes any 0s or unidentified
### species. 3
b.3 <- b[which(b$class == "migratory"), ]
b.3 <- b.3[which(b.3$Species != "Chinook"), ]
b.3 <- b.3[which(b.3$Species != "Chum"), ]
### standardize continuous variables
b.3$s.temp <- standardize(b.3$Temp.surf, centerFun = mean, scaleFun = sd)
b.3$s.sal <- standardize(b.3$Sal.surf, centerFun = mean, scaleFun = sd)
b.3$s.do <- standardize(b.3$DOmg.surf, centerFun = mean, scaleFun = sd)
b.3$s.pH <- standardize(b.3$pH.surf, centerFun = mean, scaleFun = sd)
b.3$s.J.date <- standardize(b.3$J.date, centerFun = mean, scaleFun = sd)
### Create variable j2 which is Julian day squared
b.3$j2 <- b.3$s.J.date^2
# 1 summarize by site-day
b.migra <- ddply(b.3, .(Year, J.date, s.J.date, j2, Habitat, Site, s.temp, s.sal,
s.do, s.pH), summarize, abundance = sum(abundance))
# plot abundance over julian day to see distribution:
plot(abundance ~ J.date, data = b.migra)
Habitat variables
## Add marsh habitat variables for each site
b.migra.hab <- b.migra
b.migra.hab$stnwidth <- hab2[match(b.migra.hab$Site, hab$Site), 2]
b.migra.hab$vegelev <- hab2[match(b.migra.hab$Site, hab$Site), 3]
b.migra.hab$shtdensity <- hab2[match(b.migra.hab$Site, hab$Site), 4]
b.migra.hab$shtdenhigh <- hab2[match(b.migra.hab$Site, hab$Site), 5]
b.migra.hab$tcelev <- hab2[match(b.migra.hab$Site, hab$Site), 6]
b.migra.hab$angbank <- hab2[match(b.migra.hab$Site, hab$Site), 7]
b.migra.hab$meanturb <- hab2[match(b.migra.hab$Site, hab$Site), 8]
## standardize habitat variables using 'standardize' function from robustHD
## package
b.migra.hab$stnwidth <- standardize(b.migra.hab$stnwidth, centerFun = mean,
scaleFun = sd)
b.migra.hab$stnwidth <- as.numeric(b.migra.hab$stnwidth)
b.migra.hab$vegelev <- standardize(b.migra.hab$vegelev, centerFun = mean, scaleFun = sd)
b.migra.hab$vegelev <- as.numeric(b.migra.hab$vegelev)
b.migra.hab$shtdensity <- standardize(b.migra.hab$shtdensity, centerFun = mean,
scaleFun = sd)
b.migra.hab$shtdensity <- as.numeric(b.migra.hab$shtdensity)
b.migra.hab$shtdenhigh <- standardize(b.migra.hab$shtdenhigh, centerFun = mean,
scaleFun = sd)
b.migra.hab$shtdenhigh <- as.numeric(b.migra.hab$shtdenhigh)
b.migra.hab$tcelev <- standardize(b.migra.hab$tcelev, centerFun = mean, scaleFun = sd)
b.migra.hab$tcelev <- as.numeric(b.migra.hab$tcelev)
b.migra.hab$angbank <- standardize(b.migra.hab$angbank, centerFun = mean, scaleFun = sd)
b.migra.hab$angbank <- as.numeric(b.migra.hab$angbank)
b.migra.hab$meanturb <- standardize(b.migra.hab$meanturb, centerFun = mean,
scaleFun = sd)
b.migra.hab$meanturb <- as.numeric(b.migra.hab$meanturb)VIF for collinearity of habitat variables
Assess variance inflation factors
# alias function identifies covariates that are multiples of each other - in
# this case have some habitat parameters causing issues.
alias(lm(abundance ~ s.J.date + j2 + Year + s.temp + s.sal + s.do + s.pH + stnwidth +
vegelev + shtdenhigh + shtdensity + tcelev + angbank + meanturb, data = b.migra.hab))
## Model :
## abundance ~ s.J.date + j2 + Year + s.temp + s.sal + s.do + s.pH +
## stnwidth + vegelev + shtdenhigh + shtdensity + tcelev + angbank +
## meanturb
##
## Complete :
## (Intercept) s.J.date j2
## tcelev 0 0 0
## angbank 0 0 0
## meanturb 0 0 0
## Year s.temp s.sal
## tcelev 0 0 0
## angbank 0 0 0
## meanturb 0 0 0
## s.do s.pH stnwidth
## tcelev 0 0 12768/21383
## angbank 0 0 -1190/2477
## meanturb 0 0 35115/75149
## vegelev shtdenhigh shtdensity
## tcelev 26253/53636 -207653/439893 1598375/5117629
## angbank 6410/9393 -4595/12962 -15299/268071
## meanturb -12294629/12365285 -86771/116221 -43634/67957
# We determined hab variables for Migratory should be reduced to meanturb,
# and vegelev based on significance and collinearity when doing model
# exploration.
vif(lm(abundance ~ s.J.date + j2 + Year + s.temp + s.sal + s.do + s.pH + stnwidth +
vegelev + meanturb, data = b.migra.hab)) ##all less than 5
## s.J.date j2 Year s.temp s.sal s.do s.pH
## 16.198637 6.067386 19.980166 6.905689 12.298055 7.072571 5.577524
## stnwidth vegelev meanturb
## 2.771200 1.550012 3.061611
## Pearson Corr with all vars
year <- b.migra.hab$Year
jday <- b.migra.hab$s.J.date
j2 <- b.migra.hab$j2
temp <- b.migra.hab$s.temp
do <- b.migra.hab$s.do
ph <- b.migra.hab$s.pH
sal <- b.migra.hab$s.sal
veg <- b.migra.hab$vegelev
turb <- b.migra.hab$meanturb
width <- b.migra.hab$stnwidth
shootden <- b.migra.hab$shtdensity
shtdenhigh <- b.migra.hab$shtdenhigh
elev <- b.migra.hab$tcelev
angbank <- b.migra.hab$angbank
habcovar <- cbind.data.frame(year, jday, j2, temp, do, ph, sal, veg, turb, width,
shootden, shtdenhigh, elev, angbank)
ggpairs(data = na.omit(habcovar), title = "Pearson Correlation plot for all variables in beach seine migratory fish dataset")
ggsave("/Users/Lia/Documents/Git/Fraser-salmon/Ch1.Seasonal_diversity/expl.plots/PearsonCorrAllVars_migratorymarsh.pdf")Model selection
Full model
Due to the very small sample size for other migratory fishes (# observations = 17), there were issues with the model assumptions with several parameters, and the global model was reduced to the maximum number of parameters than met assumptions and allowed the model to converge. Note that models with fewer, different parameters also converged (namely with temperature but without mean turb and stnwidth), but the AIC was worse and we felt the global model should include as many parameters as possible. The outliers in this data set caused a few parameters to always give errors: Jdate, J2 and pH. These were removed from the global model. Year had VIF over 10, so removed and AIC improved by delta 2 and VIF all were reduced to below 1.
## Migratory: full model with top habitat variables
migra <- glm.nb(abundance ~ s.do + s.sal + stnwidth + vegelev + meanturb, data = b.migra.hab,
na.action = "na.fail")
summary(migra) #AIC 104.31##
## Call:
## glm.nb(formula = abundance ~ s.do + s.sal + stnwidth + vegelev +
## meanturb, data = b.migra.hab, na.action = "na.fail", init.theta = 1.175221356,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.1737 -0.9135 -0.2309 0.2802 1.4894
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.6004 0.2689 5.952 0.00000000265 ***
## s.do 0.4460 0.3076 1.450 0.14708
## s.sal -0.8763 0.2669 -3.283 0.00103 **
## stnwidth 0.1927 0.3750 0.514 0.60726
## vegelev 0.8562 0.2958 2.895 0.00380 **
## meanturb 0.4815 0.3838 1.255 0.20965
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.1752) family taken to be 1)
##
## Null deviance: 44.223 on 16 degrees of freedom
## Residual deviance: 15.147 on 11 degrees of freedom
## AIC: 104.31
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.175
## Std. Err.: 0.431
##
## 2 x log-likelihood: -90.310vif(migra)## s.do s.sal stnwidth vegelev meanturb
## 1.353410 1.194678 1.845305 1.218116 1.743007# Check model assumptions - note that we are not actually assuming linear
# relationships with Jday, but can't tell this general function that.
sjp.glm(migra, type = "ma")## `sjp.glm()` will become deprecated in the future. Please use `plot_model()` instead.## 4 outliers removed in updated model.## # A tibble: 2 x 2
## models aic
## <chr> <dbl>
## 1 original 104.
## 2 updated 81.6
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 1.0818## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 1.1311## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 3.9819e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 1.0969## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at 1.0818## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 1.1311## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 3.9819e-17## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 1.0969## --------------------
## Check significance of terms when they entered the model...## Anova:## Warning in anova.negbin(logreg, test = "Chisq"): tests made without re-
## estimating 'theta'
## Analysis of Deviance Table
##
## Model: Negative Binomial(1.1752), link: log
##
## Response: abundance
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 16 44.223
## s.do 1 0.3035 15 43.920 0.581700
## s.sal 1 16.3559 14 27.564 0.00005249 ***
## stnwidth 1 2.6652 13 24.899 0.102563
## vegelev 1 8.3786 12 16.520 0.003797 **
## meanturb 1 1.3734 11 15.147 0.241235
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# See pseudo R2 for full model
rsquared(migra, aicc = TRUE)## Class Family Link n R.squared
## 1 negbin Negative Binomial log 17 0.657492Model selection - dredge
Then continued to dredge function to determine optimal model using AICc (migra has 17 observations).
# Generate model set
model.set.migra <- dredge(migra)
# Create model selection table
model_table.migra <- model.sel(model.set.migra)
options(scipen = 7)
names(model_table.migra) <- c("(Int)", "Mean turb", "DO", "Sal", "Width", "Veg. elev",
"df", "LL", "AICc", "delta", "weight")
kable(head(model_table.migra, n = 100), digits = 3)| (Int) | Mean turb | DO | Sal | Width | Veg. elev | df | LL | AICc | delta | weight | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 21 | 1.776 | NA | NA | -0.791 | NA | 0.752 | 4 | -47.440 | 106.212 | 0.000 | 0.29547923989 |
| 22 | 1.656 | 0.528 | NA | -0.757 | NA | 0.831 | 5 | -46.113 | 107.681 | 1.469 | 0.14177439298 |
| 5 | 2.013 | NA | NA | -0.926 | NA | NA | 3 | -50.152 | 108.150 | 1.938 | 0.11211891747 |
| 17 | 1.936 | NA | NA | NA | NA | 0.971 | 3 | -50.671 | 109.188 | 2.976 | 0.06672651142 |
| 29 | 1.700 | NA | NA | -0.760 | 0.321 | 0.721 | 5 | -46.872 | 109.198 | 2.986 | 0.06639978578 |
| 23 | 1.754 | NA | 0.267 | -0.872 | NA | 0.798 | 5 | -47.065 | 109.585 | 3.373 | 0.05471795199 |
| 18 | 1.797 | 0.595 | NA | NA | NA | 1.037 | 4 | -49.351 | 110.035 | 3.823 | 0.04369206674 |
| 13 | 1.956 | NA | NA | -0.878 | 0.407 | NA | 4 | -49.662 | 110.657 | 4.444 | 0.03202041698 |
| 6 | 1.964 | 0.398 | NA | -0.916 | NA | NA | 4 | -49.717 | 110.767 | 4.554 | 0.03031062811 |
| 24 | 1.611 | 0.611 | 0.396 | -0.873 | NA | 0.891 | 6 | -45.263 | 110.927 | 4.714 | 0.02797872063 |
| 25 | 1.868 | NA | NA | NA | 0.360 | 0.902 | 4 | -50.087 | 111.508 | 5.295 | 0.02092303537 |
| 7 | 2.004 | NA | 0.118 | -0.952 | NA | NA | 4 | -50.094 | 111.522 | 5.310 | 0.02077535264 |
| 31 | 1.641 | NA | 0.476 | -0.884 | 0.499 | 0.758 | 6 | -45.830 | 112.060 | 5.847 | 0.01587836355 |
| 30 | 1.656 | 0.530 | NA | -0.757 | -0.002 | 0.832 | 6 | -46.113 | 112.627 | 6.414 | 0.01195950591 |
| 19 | 1.938 | NA | -0.025 | NA | NA | 0.967 | 4 | -50.669 | 112.671 | 6.459 | 0.01169410704 |
| 1 | 2.360 | NA | NA | NA | NA | NA | 2 | -54.410 | 113.677 | 7.465 | 0.00707234126 |
| 15 | 1.909 | NA | 0.370 | -0.950 | 0.619 | NA | 5 | -49.208 | 113.871 | 7.658 | 0.00642054976 |
| 20 | 1.789 | 0.608 | 0.071 | NA | NA | 1.046 | 5 | -49.334 | 114.122 | 7.910 | 0.00566196129 |
| 26 | 1.797 | 0.585 | NA | NA | 0.014 | 1.034 | 5 | -49.350 | 114.156 | 7.943 | 0.00556800220 |
| 8 | 1.943 | 0.476 | 0.221 | -0.961 | NA | NA | 5 | -49.531 | 114.517 | 8.304 | 0.00464786372 |
| 14 | 1.949 | 0.219 | NA | -0.890 | 0.284 | NA | 5 | -49.574 | 114.603 | 8.391 | 0.00445190818 |
| 9 | 2.284 | NA | NA | NA | 0.607 | NA | 3 | -53.815 | 115.475 | 9.263 | 0.00287808609 |
| 27 | 1.853 | NA | 0.121 | NA | 0.392 | 0.912 | 5 | -50.043 | 115.540 | 9.328 | 0.00278593207 |
| 2 | 2.312 | 0.531 | NA | NA | NA | NA | 3 | -54.042 | 115.929 | 9.717 | 0.00229363542 |
| 3 | 2.370 | NA | -0.243 | NA | NA | NA | 3 | -54.333 | 116.513 | 10.301 | 0.00171310123 |
| 32 | 1.600 | 0.481 | 0.446 | -0.876 | 0.193 | 0.856 | 7 | -45.155 | 116.755 | 10.542 | 0.00151813805 |
| 16 | 1.905 | 0.169 | 0.355 | -0.956 | 0.514 | NA | 6 | -49.159 | 118.718 | 12.505 | 0.00056887139 |
| 10 | 2.282 | 0.128 | NA | NA | 0.539 | NA | 4 | -53.801 | 118.934 | 12.722 | 0.00051047033 |
| 11 | 2.288 | NA | -0.048 | NA | 0.593 | NA | 4 | -53.812 | 118.957 | 12.744 | 0.00050485639 |
| 28 | 1.788 | 0.583 | 0.078 | NA | 0.035 | 1.038 | 6 | -49.331 | 119.062 | 12.849 | 0.00047898391 |
| 4 | 2.321 | 0.493 | -0.143 | NA | NA | NA | 4 | -54.017 | 119.368 | 13.156 | 0.00041097803 |
| 12 | 2.287 | 0.134 | -0.061 | NA | 0.519 | NA | 5 | -53.796 | 123.046 | 16.834 | 0.00006532418 |
# Model averaging Version 1: use all models with delta AIC score less than 4
model.set.migra.4 <- get.models(model.set.migra, subset = delta < 4)
avg_model.migra.4 <- model.avg(model.set.migra.4)
summary(avg_model.migra.4)
##
## Call:
## model.avg(object = model.set.migra.4)
##
## Component model call:
## glm.nb(formula = abundance ~ <7 unique rhs>, data = b.migra.hab,
## na.action = na.fail, init.theta = <7 unique values>, link = log)
##
## Component models:
## df logLik AICc delta weight
## 35 4 -47.44 106.21 0.00 0.38
## 135 5 -46.11 107.68 1.47 0.18
## 3 3 -50.15 108.15 1.94 0.14
## 5 3 -50.67 109.19 2.98 0.09
## 345 5 -46.87 109.20 2.99 0.09
## 235 5 -47.07 109.59 3.37 0.07
## 15 4 -49.35 110.04 3.82 0.06
##
## Term codes:
## meanturb s.do s.sal stnwidth vegelev
## 1 2 3 4 5
##
## Model-averaged coefficients:
## (full average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 1.79530 0.31154 0.33729 5.323 0.0000001 ***
## s.sal -0.69562 0.38277 0.39913 1.743 0.0814 .
## vegelev 0.69379 0.40250 0.42072 1.649 0.0991 .
## meanturb 0.12919 0.27823 0.28702 0.450 0.6526
## stnwidth 0.02733 0.12363 0.12978 0.211 0.8332
## s.do 0.01872 0.10727 0.11393 0.164 0.8695
##
## (conditional average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 1.7953 0.3115 0.3373 5.323 0.0000001 ***
## s.sal -0.8102 0.2790 0.3045 2.661 0.0078 **
## vegelev 0.8101 0.3081 0.3353 2.416 0.0157 *
## meanturb 0.5439 0.3168 0.3482 1.562 0.1183
## stnwidth 0.3214 0.2919 0.3218 0.999 0.3179
## s.do 0.2672 0.3127 0.3447 0.775 0.4382
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Relative variable importance:
## s.sal vegelev meanturb stnwidth s.do
## Importance: 0.86 0.86 0.24 0.09 0.07
## N containing models: 5 6 2 1 1
migra.ci <- data.frame(confint(avg_model.migra.4, full = TRUE))
# Get pseudo R squared values for models up to delta < 4
migra.4.Rsq <- rsquared(model.set.migra.4, aicc = TRUE)
## write tables to .csv for easy comparison and plugging into appendix table
avg_model_4df.migra <- data.frame(avg_model.migra.4$msTable)
avg_model_components4.migra <- cbind(migra.4.Rsq, avg_model_4df.migra)
r = data.frame(Coeff = rownames(avg_model_4df.migra, rep(NA, length(avg_model_components4.migra))))
avg_model_components4.migra <- cbind(avg_model_components4.migra, r)
avg_model_components4.migra <- avg_model_components4.migra[, -c(6, 7)]
# write.csv(avg_model_components4.migra,
# '/Users/Lia/Documents/Git/Fraser-salmon/all.data/avg_model_components4_beachmigratory2.csv')Parameter Plot
The results of model averaging including all top ranked models up to delta AICc 4
Check fit of top model by avg – including all parameters with weight >0.5
avg.migra <- glm.nb(abundance ~ s.sal + vegelev, data = b.migra.hab, na.action = "na.fail")
summary(avg.migra) # AIC 102.88 ##
## Call:
## glm.nb(formula = abundance ~ s.sal + vegelev, data = b.migra.hab,
## na.action = "na.fail", init.theta = 0.9382988358, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.06918 -0.76237 -0.42767 0.04231 1.89455
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.7763 0.2818 6.303 0.000000000292 ***
## s.sal -0.7913 0.2645 -2.992 0.00277 **
## vegelev 0.7525 0.2910 2.586 0.00971 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.9383) family taken to be 1)
##
## Null deviance: 36.393 on 16 degrees of freedom
## Residual deviance: 16.677 on 14 degrees of freedom
## AIC: 102.88
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.938
## Std. Err.: 0.332
##
## 2 x log-likelihood: -94.879vif(avg.migra) #very low## s.sal vegelev
## 1.008413 1.008413### Plot residuals vs. fitted values
plot(fitted(avg.migra), resid(avg.migra), main = "Averaged Migratory GLM", xlab = "Fitted",
ylab = "Pearson residuals")
## q plot
qqnorm(x = b.migra.hab$abundance, y = resid(avg.migra), main = "Q-Q Averaged Migratory GLM")
qqline(resid(avg.migra), col = 2)
Resident beach seine model
response = average abundance of residents [per beach seine site]
This is 4th of 4 beach models for (1)Chinook, (2)res, (3)other migratory species and (4)resident species
Load all data - standardize variables, create subgroups from beach seine data
### Grab just resident catch - automatically removes any 0s or unidentified
### species. 4
b.4 <- b[which(b$class == "resident"), ]
### standardize continuous variables
b.4$s.temp <- standardize(b.4$Temp.surf, centerFun = mean, scaleFun = sd)
b.4$s.sal <- standardize(b.4$Sal.surf, centerFun = mean, scaleFun = sd)
b.4$s.do <- standardize(b.4$DOmg.surf, centerFun = mean, scaleFun = sd)
b.4$s.pH <- standardize(b.4$pH.surf, centerFun = mean, scaleFun = sd)
b.4$s.J.date <- standardize(b.4$J.date, centerFun = mean, scaleFun = sd)
### Create variable j2 which is Julian day squared
b.4$j2 <- b.4$s.J.date^2
# 1 summarize by site-day
b.res <- ddply(b.4, .(Year, J.date, s.J.date, j2, Habitat, Site, s.temp, s.sal,
s.do, s.pH), summarize, abundance = sum(abundance))
# plot abundance over julian day to see distribution:
plot(abundance ~ J.date, data = b.res)
Habitat variables
## Add marsh habitat variables for each site
b.res.hab <- b.res
b.res.hab$stnwidth <- hab2[match(b.res.hab$Site, hab$Site), 2]
b.res.hab$vegelev <- hab2[match(b.res.hab$Site, hab$Site), 3]
b.res.hab$shtdensity <- hab2[match(b.res.hab$Site, hab$Site), 4]
b.res.hab$shtdenhigh <- hab2[match(b.res.hab$Site, hab$Site), 5]
b.res.hab$tcelev <- hab2[match(b.res.hab$Site, hab$Site), 6]
b.res.hab$angbank <- hab2[match(b.res.hab$Site, hab$Site), 7]
b.res.hab$meanturb <- hab2[match(b.res.hab$Site, hab$Site), 8]
## standardize habitat variables using 'standardize' function from robustHD
## package
b.res.hab$stnwidth <- standardize(b.res.hab$stnwidth, centerFun = mean, scaleFun = sd)
b.res.hab$stnwidth <- as.numeric(b.res.hab$stnwidth)
b.res.hab$vegelev <- standardize(b.res.hab$vegelev, centerFun = mean, scaleFun = sd)
b.res.hab$vegelev <- as.numeric(b.res.hab$vegelev)
b.res.hab$shtdensity <- standardize(b.res.hab$shtdensity, centerFun = mean,
scaleFun = sd)
b.res.hab$shtdensity <- as.numeric(b.res.hab$shtdensity)
b.res.hab$shtdenhigh <- standardize(b.res.hab$shtdenhigh, centerFun = mean,
scaleFun = sd)
b.res.hab$shtdenhigh <- as.numeric(b.res.hab$shtdenhigh)
b.res.hab$tcelev <- standardize(b.res.hab$tcelev, centerFun = mean, scaleFun = sd)
b.res.hab$tcelev <- as.numeric(b.res.hab$tcelev)
b.res.hab$angbank <- standardize(b.res.hab$angbank, centerFun = mean, scaleFun = sd)
b.res.hab$angbank <- as.numeric(b.res.hab$angbank)
b.res.hab$meanturb <- standardize(b.res.hab$meanturb, centerFun = mean, scaleFun = sd)
b.res.hab$meanturb <- as.numeric(b.res.hab$meanturb)VIF for collinearity of habitat variables
Assess variance inflation factors
# alias function identifies covariates that are multiples of each other - in
# this case have some habitat parameters causing issues.
alias(lm(abundance ~ s.J.date + j2 + Year + s.temp + s.sal + s.do + s.pH + stnwidth +
vegelev + shtdenhigh + shtdensity + tcelev + angbank + meanturb, data = b.res.hab))
## Model :
## abundance ~ s.J.date + j2 + Year + s.temp + s.sal + s.do + s.pH +
## stnwidth + vegelev + shtdenhigh + shtdensity + tcelev + angbank +
## meanturb
##
## Complete :
## (Intercept) s.J.date j2
## tcelev 0 0 0
## angbank 0 0 0
## meanturb 0 0 0
## Year s.temp s.sal
## tcelev 0 0 0
## angbank 0 0 0
## meanturb 0 0 0
## s.do s.pH stnwidth
## tcelev 0 0 877/1622
## angbank 0 0 -293738/568293
## meanturb 0 0 8163/23578
## vegelev shtdenhigh shtdensity
## tcelev 158272651/348415924 -6656702/13163509 991/3356
## angbank 141665/188257 -948/2101 -955547/14904747
## meanturb -248443/329042 -55259/84438 -185123/372691
# We determined hab variables for Resident should be reduced to shoot
# density, angbank, and tcelev based on significance and collinearity when
# doing model exploration.
vif(lm(abundance ~ s.J.date + j2 + Year + s.temp + s.sal + s.do + s.pH + shtdensity +
tcelev + angbank, data = b.res.hab)) ##all less than 5
## s.J.date j2 Year s.temp s.sal s.do
## 3.835484 2.948591 1.521699 5.571724 2.106828 1.408997
## s.pH shtdensity tcelev angbank
## 1.238688 1.143590 1.935296 1.789858
## Pearson Corr with all vars
year <- b.res.hab$Year
jday <- b.res.hab$s.J.date
j2 <- b.res.hab$j2
temp <- b.res.hab$s.temp
do <- b.res.hab$s.do
ph <- b.res.hab$s.pH
sal <- b.res.hab$s.sal
veg <- b.res.hab$vegelev
turb <- b.res.hab$meanturb
width <- b.res.hab$stnwidth
shootden <- b.res.hab$shtdensity
shtdenhigh <- b.res.hab$shtdenhigh
elev <- b.res.hab$tcelev
angbank <- b.res.hab$angbank
habcovar <- cbind.data.frame(year, jday, j2, temp, do, ph, sal, veg, turb, width,
shootden, shtdenhigh, elev, angbank)
ggpairs(data = na.omit(habcovar), title = "Pearson Correlation plot for all variables in beach seine resident fish dataset")
ggsave("/Users/Lia/Documents/Git/Fraser-salmon/Ch1.Seasonal_diversity/expl.plots/PearsonCorrAllVars_residentmarsh.pdf")Model selection
Full model
Note shoot density came out as a stronger habitat variable than turbidity in this model (delta AIC ~5, stronger weight to top model when run without shtden or meanturb). However, it was strongly correlated to turbidity (0.779) and also highly correlated to stnwidth (0.537). To avoid removing two habitat variables, and to remain consistent with the other global models, we decided not to use shoot density. Note that temperature was very highly correlated to Julian day, J2, and salinity in this data set, and caused several of the VIF values to inflate. Once removed, all VIFs dropped to below 2 and the AIC improved, so it was removed from the global model.
## Resident: full model with top habitat variables
res <- glm.nb(abundance ~ s.J.date + j2 + Year + s.sal + s.do + s.pH + meanturb +
vegelev + stnwidth, data = b.res.hab, na.action = "na.fail")
summary(res) #749.43, VIF below 2 ##
## Call:
## glm.nb(formula = abundance ~ s.J.date + j2 + Year + s.sal + s.do +
## s.pH + meanturb + vegelev + stnwidth, data = b.res.hab, na.action = "na.fail",
## init.theta = 1.206401257, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.6124 -0.8895 -0.3893 0.0975 3.2165
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.98477 0.19314 20.631 < 2e-16 ***
## s.J.date 0.79278 0.11210 7.072 1.53e-12 ***
## j2 -0.48207 0.08914 -5.408 6.37e-08 ***
## Year 1.06309 0.24228 4.388 1.15e-05 ***
## s.sal -0.07968 0.12752 -0.625 0.532084
## s.do 0.18154 0.12509 1.451 0.146708
## s.pH -0.33039 0.11572 -2.855 0.004301 **
## meanturb 0.48323 0.13005 3.716 0.000203 ***
## vegelev 0.03291 0.13490 0.244 0.807294
## stnwidth -0.17553 0.13816 -1.270 0.203934
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.2064) family taken to be 1)
##
## Null deviance: 196.121 on 77 degrees of freedom
## Residual deviance: 85.743 on 68 degrees of freedom
## AIC: 749.43
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.206
## Std. Err.: 0.183
##
## 2 x log-likelihood: -727.426res2 <- glm.nb(abundance ~ s.J.date + j2 + Year + s.sal + s.do + s.pH + shtdensity +
vegelev, data = b.res.hab, na.action = "na.fail")
summary(res2) #744.11, VIF still below 2##
## Call:
## glm.nb(formula = abundance ~ s.J.date + j2 + Year + s.sal + s.do +
## s.pH + shtdensity + vegelev, data = b.res.hab, na.action = "na.fail",
## init.theta = 1.244095276, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.7463 -1.0202 -0.3433 0.3050 3.1930
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.96891 0.18948 20.946 < 2e-16 ***
## s.J.date 0.75583 0.11092 6.814 9.49e-12 ***
## j2 -0.47398 0.08801 -5.385 7.23e-08 ***
## Year 1.00505 0.23662 4.247 2.16e-05 ***
## s.sal -0.06401 0.12392 -0.517 0.6055
## s.do 0.22313 0.12327 1.810 0.0703 .
## s.pH -0.33959 0.11443 -2.968 0.0030 **
## shtdensity -0.48810 0.11047 -4.418 9.94e-06 ***
## vegelev -0.21295 0.11279 -1.888 0.0590 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.2441) family taken to be 1)
##
## Null deviance: 201.965 on 77 degrees of freedom
## Residual deviance: 84.981 on 69 degrees of freedom
## AIC: 744.21
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.244
## Std. Err.: 0.190
##
## 2 x log-likelihood: -724.206vif(res)## s.J.date j2 Year s.sal s.do s.pH meanturb vegelev
## 1.366545 1.911245 1.280196 1.986694 1.424802 1.269219 1.457012 1.601792
## stnwidth
## 1.675253# Check model assumptions - note that several parameters have non ideal
# relationships. Also note that there are several outliers, but these are
# real so we do not remove them (despite what this test says).
sjp.glm(res, type = "ma")## `sjp.glm()` will become deprecated in the future. Please use `plot_model()` instead.## 1 outliers removed in updated model.## # A tibble: 2 x 2
## models aic
## <chr> <dbl>
## 1 original 749.
## 2 updated 718.
## --------------------
## Check significance of terms when they entered the model...## Anova:## Warning in anova.negbin(logreg, test = "Chisq"): tests made without re-
## estimating 'theta'
## Analysis of Deviance Table
##
## Model: Negative Binomial(1.2064), link: log
##
## Response: abundance
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 77 196.121
## s.J.date 1 10.356 76 185.765 0.0012906 **
## j2 1 56.016 75 129.749 7.187e-14 ***
## Year 1 16.037 74 113.711 6.210e-05 ***
## s.sal 1 0.001 73 113.711 0.9807157
## s.do 1 1.175 72 112.536 0.2783492
## s.pH 1 10.586 71 101.950 0.0011394 **
## meanturb 1 14.465 70 87.485 0.0001428 ***
## vegelev 1 0.237 69 87.248 0.6261642
## stnwidth 1 1.505 68 85.743 0.2199672
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# See pseudo R2 for full model
rsquared(res, aicc = TRUE)## Class Family Link n R.squared
## 1 negbin Negative Binomial log 78 0.562805Model selection - dredge
Then continued to dredge function to determine optimal model using AICc (res has 17 observations).
# Generate model set
model.set.res <- dredge(res)
# Create model selection table
model_table.res <- model.sel(model.set.res)
options(scipen = 7)
names(model_table.res) <- c("(Int)", "Julian day^2", "Mean turb", "DO", "Julian day",
"pH", "Sal", "Channel width", "Veg. elev.", "Year", "df", "LL", "AICc",
"delta", "weight")
kable(head(model_table.res, n = 100), digits = 3)| (Int) | Julian day^2 | Mean turb | DO | Julian day | pH | Sal | Channel width | Veg. elev. | Year | df | LL | AICc | delta | weight | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 284 | 4.129 | -0.529 | 0.442 | NA | 0.740 | -0.301 | NA | NA | NA | 0.920 | 7 | -365.573 | 746.746 | 0.000 | 0.182920560692 |
| 288 | 4.080 | -0.520 | 0.418 | 0.168 | 0.778 | -0.349 | NA | NA | NA | 0.954 | 8 | -364.827 | 747.741 | 0.994 | 0.111266844172 |
| 348 | 4.087 | -0.529 | 0.484 | NA | 0.756 | -0.285 | NA | -0.149 | NA | 1.000 | 8 | -364.841 | 747.768 | 1.022 | 0.109746806824 |
| 352 | 4.024 | -0.516 | 0.461 | 0.188 | 0.798 | -0.335 | NA | -0.165 | NA | 1.046 | 9 | -363.912 | 748.471 | 1.725 | 0.077230890598 |
| 316 | 4.074 | -0.491 | 0.452 | NA | 0.736 | -0.301 | -0.087 | NA | NA | 0.959 | 8 | -365.339 | 748.765 | 2.019 | 0.066654400522 |
| 412 | 4.118 | -0.530 | 0.430 | NA | 0.742 | -0.303 | NA | NA | -0.043 | 0.945 | 8 | -365.501 | 749.089 | 2.343 | 0.056694729815 |
| 320 | 4.023 | -0.480 | 0.429 | 0.170 | 0.774 | -0.349 | -0.090 | NA | NA | 0.991 | 9 | -364.575 | 749.796 | 3.050 | 0.039808074651 |
| 416 | 4.055 | -0.519 | 0.394 | 0.186 | 0.784 | -0.357 | NA | NA | -0.074 | 0.997 | 9 | -364.621 | 749.889 | 3.143 | 0.038007892528 |
| 380 | 4.043 | -0.497 | 0.491 | NA | 0.751 | -0.286 | -0.075 | -0.144 | NA | 1.029 | 9 | -364.659 | 749.966 | 3.219 | 0.036577090674 |
| 476 | 4.091 | -0.528 | 0.502 | NA | 0.757 | -0.280 | NA | -0.172 | 0.042 | 0.990 | 9 | -364.793 | 750.233 | 3.487 | 0.031996820852 |
| 268 | 4.085 | -0.488 | 0.465 | NA | 0.734 | NA | NA | NA | NA | 1.027 | 6 | -368.526 | 750.235 | 3.489 | 0.031968121608 |
| 332 | 4.021 | -0.488 | 0.517 | NA | 0.751 | NA | NA | -0.198 | NA | 1.142 | 7 | -367.362 | 750.324 | 3.578 | 0.030577406544 |
| 384 | 3.984 | -0.485 | 0.468 | 0.187 | 0.793 | -0.335 | -0.073 | -0.159 | NA | 1.070 | 10 | -363.740 | 750.763 | 4.017 | 0.024549895481 |
| 480 | 4.026 | -0.516 | 0.466 | 0.186 | 0.798 | -0.333 | NA | -0.172 | 0.012 | 1.043 | 10 | -363.908 | 751.100 | 4.354 | 0.020740718087 |
| 444 | 4.072 | -0.495 | 0.444 | NA | 0.737 | -0.303 | -0.080 | NA | -0.027 | 0.972 | 9 | -365.313 | 751.273 | 4.527 | 0.019020837136 |
| 448 | 4.014 | -0.486 | 0.408 | 0.183 | 0.778 | -0.356 | -0.074 | NA | -0.058 | 1.019 | 10 | -364.459 | 752.201 | 5.454 | 0.011964306605 |
| 336 | 3.976 | -0.475 | 0.509 | 0.104 | 0.774 | NA | NA | -0.214 | NA | 1.186 | 8 | -367.096 | 752.278 | 5.532 | 0.011509420075 |
| 300 | 4.030 | -0.451 | 0.477 | NA | 0.732 | NA | -0.080 | NA | NA | 1.068 | 7 | -368.347 | 752.294 | 5.548 | 0.011415914866 |
| 508 | 4.042 | -0.490 | 0.520 | NA | 0.752 | -0.278 | -0.089 | -0.179 | 0.065 | 1.017 | 10 | -364.554 | 752.391 | 5.645 | 0.010877766864 |
| 272 | 4.063 | -0.481 | 0.458 | 0.064 | 0.748 | NA | NA | NA | NA | 1.048 | 7 | -368.425 | 752.449 | 5.703 | 0.010564950833 |
| 460 | 4.033 | -0.488 | 0.549 | NA | 0.755 | NA | NA | -0.239 | 0.080 | 1.111 | 8 | -367.202 | 752.490 | 5.744 | 0.010350615802 |
| 364 | 3.980 | -0.459 | 0.525 | NA | 0.748 | NA | -0.065 | -0.195 | NA | 1.172 | 8 | -367.233 | 752.553 | 5.806 | 0.010033063566 |
| 396 | 4.073 | -0.488 | 0.456 | NA | 0.734 | NA | NA | NA | -0.034 | 1.052 | 7 | -368.485 | 752.571 | 5.824 | 0.009942402590 |
| 512 | 3.985 | -0.482 | 0.483 | 0.182 | 0.793 | -0.330 | -0.080 | -0.176 | 0.033 | 1.063 | 11 | -363.713 | 753.426 | 6.680 | 0.006482096312 |
| 492 | 3.980 | -0.445 | 0.572 | NA | 0.753 | NA | -0.094 | -0.248 | 0.111 | 1.138 | 9 | -366.955 | 754.557 | 7.810 | 0.003683410824 |
| 304 | 4.007 | -0.443 | 0.469 | 0.064 | 0.746 | NA | -0.081 | NA | NA | 1.088 | 8 | -368.243 | 754.572 | 7.826 | 0.003654809769 |
| 368 | 3.938 | -0.447 | 0.517 | 0.102 | 0.771 | NA | -0.063 | -0.210 | NA | 1.212 | 9 | -366.976 | 754.599 | 7.853 | 0.003606288205 |
| 464 | 3.991 | -0.476 | 0.536 | 0.094 | 0.775 | NA | NA | -0.246 | 0.066 | 1.157 | 9 | -366.992 | 754.630 | 7.884 | 0.003550787903 |
| 428 | 4.029 | -0.453 | 0.472 | NA | 0.732 | NA | -0.075 | NA | -0.015 | 1.077 | 8 | -368.340 | 754.768 | 8.021 | 0.003314812243 |
| 400 | 4.040 | -0.479 | 0.443 | 0.077 | 0.750 | NA | NA | NA | -0.050 | 1.088 | 8 | -368.345 | 754.777 | 8.030 | 0.003299759334 |
| 496 | 3.945 | -0.437 | 0.559 | 0.087 | 0.772 | NA | -0.089 | -0.254 | 0.096 | 1.177 | 10 | -366.775 | 756.834 | 10.087 | 0.001179934982 |
| 432 | 4.001 | -0.447 | 0.459 | 0.072 | 0.747 | NA | -0.071 | NA | -0.030 | 1.108 | 9 | -368.217 | 757.080 | 10.334 | 0.001042883619 |
| 414 | 4.141 | -0.542 | NA | 0.294 | 0.761 | -0.419 | NA | NA | -0.193 | 0.984 | 8 | -369.676 | 757.439 | 10.693 | 0.000871613823 |
| 286 | 4.235 | -0.554 | NA | 0.262 | 0.744 | -0.412 | NA | NA | NA | 0.856 | 7 | -371.167 | 757.934 | 11.188 | 0.000680572339 |
| 282 | 4.317 | -0.572 | NA | NA | 0.696 | -0.338 | NA | NA | NA | 0.820 | 6 | -372.714 | 758.611 | 11.864 | 0.000485219072 |
| 410 | 4.250 | -0.568 | NA | NA | 0.709 | -0.343 | NA | NA | -0.163 | 0.928 | 7 | -371.651 | 758.902 | 12.155 | 0.000419546007 |
| 28 | 4.642 | -0.555 | 0.396 | NA | 0.631 | -0.371 | NA | NA | NA | NA | 6 | -372.956 | 759.094 | 12.348 | 0.000381018984 |
| 478 | 4.144 | -0.542 | NA | 0.289 | 0.759 | -0.423 | NA | 0.043 | -0.209 | 0.973 | 9 | -369.623 | 759.893 | 13.147 | 0.000255537591 |
| 446 | 4.145 | -0.546 | NA | 0.294 | 0.762 | -0.419 | 0.008 | NA | -0.194 | 0.982 | 9 | -369.674 | 759.996 | 13.249 | 0.000242779941 |
| 350 | 4.221 | -0.553 | NA | 0.272 | 0.750 | -0.409 | NA | -0.050 | NA | 0.883 | 8 | -371.074 | 760.235 | 13.489 | 0.000215355967 |
| 318 | 4.217 | -0.539 | NA | 0.263 | 0.744 | -0.413 | -0.033 | NA | NA | 0.866 | 8 | -371.138 | 760.362 | 13.616 | 0.000202111334 |
| 314 | 4.300 | -0.559 | NA | NA | 0.696 | -0.338 | -0.030 | NA | NA | 0.832 | 7 | -372.689 | 760.978 | 14.232 | 0.000148554064 |
| 32 | 4.622 | -0.549 | 0.373 | 0.097 | 0.652 | -0.404 | NA | NA | NA | NA | 7 | -372.695 | 760.991 | 14.245 | 0.000147604798 |
| 346 | 4.313 | -0.573 | NA | NA | 0.699 | -0.336 | NA | -0.022 | NA | 0.832 | 7 | -372.696 | 760.993 | 14.246 | 0.000147479176 |
| 474 | 4.252 | -0.566 | NA | NA | 0.706 | -0.349 | NA | 0.063 | -0.186 | 0.912 | 8 | -371.537 | 761.162 | 14.415 | 0.000135535662 |
| 156 | 4.637 | -0.552 | 0.415 | NA | 0.634 | -0.357 | NA | NA | 0.066 | NA | 7 | -372.809 | 761.217 | 14.471 | 0.000131818091 |
| 442 | 4.250 | -0.568 | NA | NA | 0.709 | -0.343 | -0.001 | NA | -0.163 | 0.929 | 8 | -371.651 | 761.389 | 14.642 | 0.000120992507 |
| 60 | 4.654 | -0.568 | 0.394 | NA | 0.635 | -0.370 | 0.028 | NA | NA | NA | 7 | -372.936 | 761.472 | 14.725 | 0.000116061701 |
| 92 | 4.642 | -0.555 | 0.401 | NA | 0.632 | -0.372 | NA | -0.016 | NA | NA | 7 | -372.948 | 761.496 | 14.749 | 0.000114689361 |
| 266 | 4.314 | -0.536 | NA | NA | 0.680 | NA | NA | NA | NA | 0.869 | 5 | -375.831 | 762.495 | 15.748 | 0.000069592784 |
| 510 | 4.147 | -0.545 | NA | 0.289 | 0.759 | -0.423 | 0.006 | 0.042 | -0.209 | 0.971 | 10 | -369.622 | 762.528 | 15.781 | 0.000068455690 |
| 382 | 4.207 | -0.542 | NA | 0.272 | 0.749 | -0.410 | -0.027 | -0.047 | NA | 0.890 | 9 | -371.055 | 762.756 | 16.010 | 0.000061064849 |
| 394 | 4.226 | -0.531 | NA | NA | 0.685 | NA | NA | NA | -0.161 | 1.016 | 6 | -374.899 | 762.980 | 16.234 | 0.000054592919 |
| 160 | 4.622 | -0.548 | 0.391 | 0.086 | 0.652 | -0.390 | NA | NA | 0.050 | NA | 8 | -372.615 | 763.316 | 16.570 | 0.000046148090 |
| 378 | 4.297 | -0.560 | NA | NA | 0.698 | -0.337 | -0.029 | -0.020 | NA | 0.842 | 8 | -372.674 | 763.435 | 16.688 | 0.000043493776 |
| 96 | 4.622 | -0.550 | 0.379 | 0.099 | 0.654 | -0.406 | NA | -0.024 | NA | NA | 8 | -372.678 | 763.444 | 16.697 | 0.000043298394 |
| 220 | 4.638 | -0.553 | 0.447 | NA | 0.640 | -0.353 | NA | -0.076 | 0.108 | NA | 8 | -372.682 | 763.452 | 16.705 | 0.000043123246 |
| 64 | 4.632 | -0.560 | 0.372 | 0.096 | 0.656 | -0.402 | 0.022 | NA | NA | NA | 8 | -372.683 | 763.452 | 16.706 | 0.000043112910 |
| 188 | 4.640 | -0.555 | 0.414 | NA | 0.635 | -0.357 | 0.006 | NA | 0.065 | NA | 8 | -372.808 | 763.703 | 16.956 | 0.000038043331 |
| 506 | 4.252 | -0.566 | NA | NA | 0.706 | -0.349 | -0.001 | 0.063 | -0.186 | 0.912 | 9 | -371.537 | 763.722 | 16.975 | 0.000037683780 |
| 124 | 4.655 | -0.569 | 0.399 | NA | 0.637 | -0.370 | 0.029 | -0.018 | NA | NA | 8 | -372.926 | 763.939 | 17.193 | 0.000033795603 |
| 270 | 4.264 | -0.520 | NA | 0.145 | 0.705 | NA | NA | NA | NA | 0.901 | 6 | -375.383 | 763.948 | 17.202 | 0.000033644199 |
| 398 | 4.142 | -0.506 | NA | 0.187 | 0.715 | NA | NA | NA | -0.187 | 1.083 | 7 | -374.176 | 763.952 | 17.206 | 0.000033580163 |
| 330 | 4.298 | -0.538 | NA | NA | 0.685 | NA | NA | -0.059 | NA | 0.907 | 6 | -375.716 | 764.614 | 17.868 | 0.000024117160 |
| 298 | 4.312 | -0.534 | NA | NA | 0.680 | NA | -0.004 | NA | NA | 0.871 | 6 | -375.830 | 764.843 | 18.097 | 0.000021504137 |
| 426 | 4.246 | -0.547 | NA | NA | 0.686 | NA | 0.036 | NA | -0.167 | 1.003 | 7 | -374.865 | 765.330 | 18.584 | 0.000016860157 |
| 458 | 4.226 | -0.530 | NA | NA | 0.684 | NA | NA | 0.009 | -0.164 | 1.014 | 7 | -374.896 | 765.393 | 18.646 | 0.000016340838 |
| 224 | 4.622 | -0.549 | 0.422 | 0.084 | 0.658 | -0.384 | NA | -0.073 | 0.090 | NA | 9 | -372.499 | 765.645 | 18.898 | 0.000014406394 |
| 12 | 4.694 | -0.516 | 0.382 | NA | 0.658 | NA | NA | NA | NA | NA | 5 | -377.426 | 765.686 | 18.940 | 0.000014109604 |
| 192 | 4.624 | -0.551 | 0.390 | 0.086 | 0.653 | -0.390 | 0.006 | NA | 0.049 | NA | 9 | -372.614 | 765.875 | 19.128 | 0.000012841836 |
| 334 | 4.231 | -0.518 | NA | 0.167 | 0.714 | NA | NA | -0.085 | NA | 0.960 | 7 | -375.152 | 765.905 | 19.158 | 0.000012650753 |
| 128 | 4.633 | -0.561 | 0.378 | 0.098 | 0.658 | -0.404 | 0.024 | -0.025 | NA | NA | 9 | -372.663 | 765.974 | 19.227 | 0.000012221557 |
| 252 | 4.637 | -0.552 | 0.448 | NA | 0.640 | -0.353 | -0.003 | -0.077 | 0.108 | NA | 9 | -372.682 | 766.012 | 19.265 | 0.000011991232 |
| 430 | 4.167 | -0.529 | NA | 0.191 | 0.717 | NA | 0.051 | NA | -0.197 | 1.069 | 8 | -374.112 | 766.312 | 19.565 | 0.000010321217 |
| 302 | 4.263 | -0.519 | NA | 0.145 | 0.705 | NA | -0.002 | NA | NA | 0.902 | 7 | -375.382 | 766.365 | 19.619 | 0.000010049514 |
| 462 | 4.140 | -0.506 | NA | 0.189 | 0.716 | NA | NA | -0.013 | -0.182 | 1.087 | 8 | -374.171 | 766.430 | 19.683 | 0.000009729040 |
| 140 | 4.680 | -0.512 | 0.426 | NA | 0.668 | NA | NA | NA | 0.143 | NA | 6 | -376.761 | 766.705 | 19.959 | 0.000008477039 |
| 26 | 4.767 | -0.593 | NA | NA | 0.596 | -0.366 | NA | NA | NA | NA | 5 | -377.950 | 766.734 | 19.987 | 0.000008357446 |
| 362 | 4.297 | -0.537 | NA | NA | 0.685 | NA | -0.001 | -0.059 | NA | 0.908 | 7 | -375.715 | 767.031 | 20.285 | 0.000007203186 |
| 30 | 4.712 | -0.576 | NA | 0.199 | 0.633 | -0.433 | NA | NA | NA | NA | 6 | -376.969 | 767.121 | 20.375 | 0.000006884766 |
| 44 | 4.733 | -0.558 | 0.374 | NA | 0.668 | NA | 0.084 | NA | NA | NA | 6 | -377.256 | 767.695 | 20.948 | 0.000005168447 |
| 490 | 4.247 | -0.547 | NA | NA | 0.685 | NA | 0.037 | 0.011 | -0.171 | 1.000 | 8 | -374.862 | 767.811 | 21.065 | 0.000004877129 |
| 16 | 4.702 | -0.521 | 0.391 | -0.038 | 0.649 | NA | NA | NA | NA | NA | 6 | -377.388 | 767.959 | 21.213 | 0.000004528053 |
| 76 | 4.694 | -0.516 | 0.383 | NA | 0.658 | NA | NA | -0.005 | NA | NA | 6 | -377.426 | 768.035 | 21.288 | 0.000004360744 |
| 256 | 4.622 | -0.548 | 0.422 | 0.084 | 0.658 | -0.384 | -0.001 | -0.073 | 0.091 | NA | 10 | -372.499 | 768.281 | 21.535 | 0.000003855348 |
| 366 | 4.234 | -0.520 | NA | 0.167 | 0.714 | NA | 0.004 | -0.085 | NA | 0.958 | 8 | -375.152 | 768.391 | 21.644 | 0.000003650020 |
| 204 | 4.679 | -0.513 | 0.478 | NA | 0.679 | NA | NA | -0.127 | 0.210 | NA | 7 | -376.442 | 768.485 | 21.738 | 0.000003482261 |
| 90 | 4.757 | -0.588 | NA | NA | 0.593 | -0.364 | NA | 0.072 | NA | NA | 6 | -377.773 | 768.728 | 21.982 | 0.000003082673 |
| 154 | 4.765 | -0.594 | NA | NA | 0.595 | -0.378 | NA | NA | -0.055 | NA | 6 | -377.839 | 768.861 | 22.115 | 0.000002885089 |
| 494 | 4.165 | -0.528 | NA | 0.193 | 0.718 | NA | 0.050 | -0.013 | -0.192 | 1.073 | 9 | -374.108 | 768.863 | 22.117 | 0.000002881381 |
| 58 | 4.794 | -0.622 | NA | NA | 0.603 | -0.362 | 0.062 | NA | NA | NA | 6 | -377.859 | 768.901 | 22.155 | 0.000002827811 |
| 144 | 4.690 | -0.518 | 0.444 | -0.063 | 0.652 | NA | NA | NA | 0.150 | NA | 7 | -376.659 | 768.918 | 22.172 | 0.000002803762 |
| 158 | 4.703 | -0.575 | NA | 0.209 | 0.633 | -0.452 | NA | NA | -0.078 | NA | 7 | -376.750 | 769.100 | 22.354 | 0.000002559527 |
| 172 | 4.689 | -0.521 | 0.423 | NA | 0.669 | NA | 0.018 | NA | 0.136 | NA | 7 | -376.755 | 769.110 | 22.364 | 0.000002547133 |
| 94 | 4.709 | -0.575 | NA | 0.190 | 0.629 | -0.429 | NA | 0.049 | NA | NA | 7 | -376.885 | 769.370 | 22.624 | 0.000002236313 |
| 62 | 4.734 | -0.600 | NA | 0.196 | 0.638 | -0.429 | 0.050 | NA | NA | NA | 7 | -376.911 | 769.421 | 22.675 | 0.000002180297 |
| 299 | 3.558 | NA | 0.599 | NA | 0.669 | NA | -0.439 | NA | NA | 1.192 | 6 | -378.375 | 769.932 | 23.186 | 0.000001688567 |
| 48 | 4.740 | -0.563 | 0.383 | -0.039 | 0.658 | NA | 0.084 | NA | NA | NA | 7 | -377.217 | 770.033 | 23.287 | 0.000001605295 |
| 108 | 4.734 | -0.559 | 0.377 | NA | 0.668 | NA | 0.086 | -0.014 | NA | NA | 7 | -377.251 | 770.102 | 23.355 | 0.000001551472 |
| 80 | 4.702 | -0.521 | 0.392 | -0.038 | 0.649 | NA | NA | -0.002 | NA | NA | 7 | -377.388 | 770.376 | 23.630 | 0.000001352465 |
# Model averaging Version 1: use all models with delta AIC score less than 4
model.set.res.4 <- get.models(model.set.res, subset = delta < 4)
avg_model.res.4 <- model.avg(model.set.res.4)
summary(avg_model.res.4)
##
## Call:
## model.avg(object = model.set.res.4)
##
## Component model call:
## glm.nb(formula = abundance ~ <12 unique rhs>, data = b.res.hab,
## na.action = na.fail, init.theta = <12 unique values>, link = log)
##
## Component models:
## df logLik AICc delta weight
## 12459 7 -365.57 746.75 0.00 0.22
## 123459 8 -364.83 747.74 0.99 0.14
## 124579 8 -364.84 747.77 1.02 0.13
## 1234579 9 -363.91 748.47 1.72 0.09
## 124569 8 -365.34 748.77 2.02 0.08
## 124589 8 -365.50 749.09 2.34 0.07
## 1234569 9 -364.57 749.80 3.05 0.05
## 1234589 9 -364.62 749.89 3.14 0.05
## 1245679 9 -364.66 749.97 3.22 0.04
## 1245789 9 -364.79 750.23 3.49 0.04
## 1249 6 -368.53 750.23 3.49 0.04
## 12479 7 -367.36 750.32 3.58 0.04
##
## Term codes:
## j2 meanturb s.do s.J.date s.pH s.sal stnwidth vegelev
## 1 2 3 4 5 6 7 8
## Year
## 9
##
## Model-averaged coefficients:
## (full average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 4.081422 0.180832 0.183840 22.201 < 2e-16 ***
## j2 -0.516061 0.073616 0.074828 6.897 < 2e-16 ***
## meanturb 0.451576 0.117970 0.119888 3.767 0.000165 ***
## s.J.date 0.757814 0.107372 0.109167 6.942 < 2e-16 ***
## s.pH -0.288963 0.138021 0.139470 2.072 0.038277 *
## Year 0.979812 0.237464 0.241384 4.059 0.0000493 ***
## s.do 0.057912 0.109888 0.110713 0.523 0.600919
## stnwidth -0.056596 0.105597 0.106453 0.532 0.594969
## s.sal -0.014969 0.061265 0.062036 0.241 0.809324
## vegelev -0.004848 0.050826 0.051565 0.094 0.925098
##
## (conditional average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 4.08142 0.18083 0.18384 22.201 < 2e-16 ***
## j2 -0.51606 0.07362 0.07483 6.897 < 2e-16 ***
## meanturb 0.45158 0.11797 0.11989 3.767 0.000165 ***
## s.J.date 0.75781 0.10737 0.10917 6.942 < 2e-16 ***
## s.pH -0.31303 0.11446 0.11635 2.690 0.007137 **
## Year 0.97981 0.23746 0.24138 4.059 0.0000493 ***
## s.do 0.17689 0.12585 0.12804 1.382 0.167114
## stnwidth -0.16090 0.12214 0.12424 1.295 0.195287
## s.sal -0.08513 0.12399 0.12615 0.675 0.499785
## vegelev -0.03112 0.12557 0.12749 0.244 0.807126
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Relative variable importance:
## j2 meanturb s.J.date Year s.pH stnwidth s.do s.sal
## Importance: 1.00 1.00 1.00 1.00 0.92 0.35 0.33 0.18
## N containing models: 12 12 12 12 10 5 4 3
## vegelev
## Importance: 0.16
## N containing models: 3
res.ci <- data.frame(confint(avg_model.res.4, full = TRUE))
# Get pseudo R squared values for models up to delta < 4
res.4.Rsq <- rsquared(model.set.res.4, aicc = TRUE)
## write tables to .csv for easy comparison and plugging into appendix table
avg_model_4df.res <- data.frame(avg_model.res.4$msTable)
avg_model_components4.res <- cbind(res.4.Rsq, avg_model_4df.res)
r = data.frame(Coeff = rownames(avg_model_4df.res, rep(NA, length(avg_model_components4.res))))
avg_model_components4.res <- cbind(avg_model_components4.res, r)
avg_model_components4.res <- avg_model_components4.res[, -c(6, 7)]
# write.csv(avg_model_components4.res,
# '/Users/Lia/Documents/Git/Fraser-salmon/all.data/avg_model_components4_beachresident.csv')Parameter Plot
The results of model averaging including all top ranked models up to delta AICc 4
Check fit of top model by avg – including all parameters with weight >0.5
avg.res <- glm.nb(abundance ~ Year + s.J.date + j2 + meanturb + s.pH, data = b.res.hab,
na.action = "na.fail")
summary(avg.res)##
## Call:
## glm.nb(formula = abundance ~ Year + s.J.date + j2 + meanturb +
## s.pH, data = b.res.hab, na.action = "na.fail", init.theta = 1.152490818,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.6782 -0.9139 -0.3890 0.1030 3.2421
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.12855 0.17335 23.816 < 2e-16 ***
## Year 0.92038 0.22756 4.045 5.24e-05 ***
## s.J.date 0.74024 0.10212 7.249 4.21e-13 ***
## j2 -0.52937 0.06791 -7.795 6.43e-15 ***
## meanturb 0.44222 0.11029 4.010 6.08e-05 ***
## s.pH -0.30052 0.10907 -2.755 0.00586 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.1525) family taken to be 1)
##
## Null deviance: 187.737 on 77 degrees of freedom
## Residual deviance: 85.784 on 72 degrees of freedom
## AIC: 745.15
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.152
## Std. Err.: 0.174
##
## 2 x log-likelihood: -731.146vif(avg.res)## Year s.J.date j2 meanturb s.pH
## 1.086656 1.103259 1.089915 1.007426 1.082319### Plot residuals vs. fitted values
plot(fitted(avg.res), resid(avg.res), main = "Averaged Resident GLM", xlab = "Fitted",
ylab = "Pearson residuals")
## q plot
qqnorm(x = b.res.hab$abundance, y = resid(avg.res), main = "Q-Q Averaged Resident GLM")
qqline(resid(avg.res), col = 2)